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Professor  Mikhail  D.  Lukin  Michael  John  Gullans 

Controlling  Atomic,  Solid-State  and  Hybrid  Systems  for 
Quantum  Information  Processing 

Abstract 

Quantum  information  science  involves  the  use  of  precise  control  over  quantum 
systems  to  explore  new  technologies.  However,  as  quantum  systems  are  scaled  up 
they  require  an  ever  deeper  understanding  of  many-body  physics  to  achieve  the  re¬ 
quired  degree  of  control.  Current  experiments  are  entering  a  regime  which  requires 
active  control  of  a  mesoscopic  number  of  coupled  quantum  systems  or  quantum  bits 
(qubits).  This  thesis  describes  several  approaches  to  this  goal  and  shows  how  meso¬ 
scopic  quantum  systems  can  be  controlled  and  utilized  for  quantum  information  tasks. 

The  first  system  we  consider  is  the  nuclear  spin  environment  of  GaAs  double  quan¬ 
tum  dots  containing  two  electrons.  We  show  that  the  through  appropriate  control  of 
dynamic  nuclear  polarization  one  can  prepare  the  nuclear  spin  environment  in  three 
distinct  collective  quantum  states  which  are  useful  for  quantum  information  process¬ 
ing  with  electron  spin  qubits.  We  then  investigate  a  hybrid  system  in  which  an  optical 
lattice  is  formed  in  the  near  field  scattering  off  an  array  of  metallic  nanoparticles  by 
utilizing  the  plasmonic  resonance  of  the  nanoparticles.  We  show  that  such  a  system 
would  realize  new  regimes  of  dense,  ultra-cold  quantum  matter  and  can  be  used  to 
create  a  quantum  network  of  atoms  and  plasmons.  Finally  we  investigate  quantum 
nonlinear  optical  systems.  We  show  that  the  intrinsic  nonlinearity  for  plasmons  in 
graphene  can  be  large  enough  to  make  a  quantum  gate  for  single  photons.  We  also 

iii 


Abstract 


consider  two  nonlinear  optical  systems  based  on  ultracold  gases  of  atoms.  In  one 
case,  we  theoretically  analyze  an  all-optical  single  photon  switch  using  cavity  quan¬ 
tum  electrodynamics  (QED)  and  slow  light.  In  the  second  case,  we  study  few  photon 
physics  in  strongly  interacting  Rydberg  polariton  systems,  where  we  demonstrate  the 
existence  of  two  and  three  photon  bound  states  and  study  their  properties. 
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Chapter  1 
Introduction 

1.1  Motivation 

Since  its  inception  in  the  80’s  and  90’s  quantum  information  science  has  developed 
into  a  mature  held  whose  central  goal  is  to  develop  new  technologies  based  on  the  pre¬ 
cise  control  of  quantum  systems.  Realization  of  this  goal  requires  contributions  from 
many  fields  of  science  and  engineering  including  physics,  materials  science,  computer 
science,  chemistry  and,  even,  biology.  Broadly  speaking  the  applications  for  such 
quantum  systems  fall  into  two  categories:  information  science,  i.e.  computation  and 
communication,  and  measurement  science,  i.e.  improved  (broadly  defined)  sensors 
and  precision.  On  the  surface  these  two  sets  of  applications  seem  unrelated,  however, 
they  are  intricately  linked  in  quantum  science.  As  quantum  information  systems  are 
pushed  to  their  limits  in  terms  of  complexity  they  require  increasing  precision  to  char¬ 
acterize  and  operate.  In  addition,  quantum  systems  developed  for  information  science 
are  so  well  isolated  and  controllable  that  developing  them  into  precision  sensors  is  a 
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natural  application.  At  the  same  time,  quantum  systems  developed  for  their  metro¬ 
logical  applications  have  become  good  candidates  for  the  building  blocks  of  quantum 
computers. 

This  thesis  will  focus  on  the  use  of  quantum  systems  for  applications  in  informa¬ 
tion  science.  The  challenges  in  this  held  include,  first,  scaling  up  the  quantum  systems 
and,  then,  achieving  sufficient  control  to  utilize  them  for  information  science.  This  is 
difficult  on  a  technological  level  as  it  requires  the  development  of  new  systems  with 
improved  control,  but  also  on  a  fundamental  level  because  the  use  of  such  systems 
requires  a  deep  understanding  of  the  many-body  physics  of  interacting  quantum  sys¬ 
tems.  Gaining  such  understanding  is,  perhaps,  the  primary  goal  in  theoretical  efforts 
for  quantum  information  science.  A  task  which  is  often  complicated  by  the  fact  that 
the  systems  under  consideration  are  fundamentally  out  of  equilibrium.  The  interplay 
between  non-equilibrium  many-body  physics  and  quantum  information  science  is  a 
central  theme  in  this  thesis.  In  what  follows  we  explore  a  range  of  physical  systems 
currently  being  pursued  for  quantum  information  applications  with  the  goal  of  har¬ 
nessing  their  many-body  behavior  to  achieve  new  applications  as  well  as  a  deeper 
understanding  of  quantum  physics. 

1.2  Mesoscopic  Quantum  Systems 

1.2.1  Electron  and  Nuclear  Spins  in  Solids 

A  promising  candidate  for  a  qubit,  the  fundamental  building  block  of  a  quantum 
computer,  is  the  spin  of  an  electron.  The  electron  can  be  bound  to  an  atom  or  ion  in 
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free  space  or  confined  in  a  solid  state  environment.  In  the  former  case  the  electron 
is  well  isolated  with  long  coherence  times;  however,  the  fabrication  and  manipulation 
of  such  systems  is  cumbersome  making  it  difficult  to  scale  them  up  to  many  qubits. 
On  the  other  hand,  trapping  electrons  in  solids  holds  promising  potential  for  scaling 
up  to  a  full  size  quantum  computer,  but  has  the  tradeoff  that  the  electrons  interacts 
strongly  with  their  host  environment.  Nevertheless  there  are  several  condensed  matter 
systems  where  single  electrons  can  be  well  enough  isolated  from  their  environment 
that  they  have  coherence  properties  comparable  to  single  atoms  or  ions. 

Solid  state  spin  qubits  generally  arise  from  the  electron  spin  of  an  impurity  atom 
or  a  quantum  dot.  Notable  examples  of  impurity  systems  are  Nitrogen  Vacancy 
(NV)  centers  in  diamond  (Jelezko  and  Wrachtrup,  2006)  and  phosphorous  donors 
in  silicon  (Zwanenburg  et  al.,  2013).  Such  impurities,  to  a  large  degree,  behaves 
as  single  trapped  atoms.  Quantum  dots  are  artificially  trapped  electrons  which  are 
confined  by  a  material  interface,  sometimes  in  combination  with  electric  gates  (Hanson 
et  al.,  2007).  Many  properties  of  quantum  dot  systems  can  also  be  explained  by 
treating  them  as  atoms,  but  the  typical  confinement  energies  (from  meV  to  eV)  and 
length  scales  (from  several  microns  to  a  few  nanometers)  vary  more  dramatically  than 
impurity  based  qubits.  In  this  thesis  we  focus  on  electrically  gated  double  quantum 
dots  in  type  III-V  semiconductors  whose  properties  are  illustrated  in  Fig.  1.1. 

Many  fundamental  quantum  operations  have  been  demonstrated  for  such  double 
quantum  dot  systems  including,  initialization,  readout  and  single  qubit  operations 
(Petta  et  ah,  2005)  and  two  qubit  entanglement  (Shulman  et  ah,  2012).  However,  a 
ubiquitous  problem  with  these  systems  is  that  the  electron  spins  interact  with  the 
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Figure  1.1:  a)  Energy  diagram  of  a  GaAs  double  quantum  dot,  where  (n,  m )  refers 
to  the  charge  occupancy  of  the  two  dots  and  S(Tm)  refers  to  the  two  electron  spin 
state  in  a  singlet  (triplet  m)  state.  J  is  the  exchange  splitting  between  the  two  singlet 
states,  Bext  is  the  magnetic  field  and  e  is  the  voltage  difference  between  the  left  and 
right  dots.  The  qubit  states  are  formed  from  the  (1,1)5  and  (1,1)T0  state,  b)  SEM 
image  of  a  double  quantum  dot  in  GaAs  (adapted  from  (Petta  et  al,  2005)).  c)  A 
double  quantum  dot  with  two  electrons  interacting  with  a  large  number  of  lattice 
nuclear  spins. 
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nuclear  spin  of  the  GaAs  host  lattice  through  the  hyperhne  interaction.  As  the  scale 
of  these  quantum  dots  is  several  hundred  nanometers,  the  electrons  interact  with  on 
the  order  of  105  —  106  nuclei  as  illustrated  in  Fig.  l.lbc.  A  surprising  feature  of 
the  nuclear  spin  environment  is  that,  due  to  long  nuclear  spin  coherence  times,  one 
can  use  the  nuclear  spins  as  a  resource  for  quantum  control  of  the  double  dot  qubit 
(Foletti  et  ah,  2009).  In  this  thesis  we  explore  theoretically  how  to  achieve  this  control 
through  dynamic  nuclear  polarization  of  the  nuclear  spins. 

1.2.2  Optical  Lattices  for  Ultracold  Atoms 

Building  a  general  purpose  quantum  computer  remains  an  outstanding  challenge. 
A  more  immediate  goal  is  to  build  a  quantum  simulator,  which  is  a  device  that  can 
solve  the  quantum  dynamics  of  an  interacting,  many-body  Hamiltonian.  A  powerful 
realization  of  a  quantum  simulator  is  an  ensemble  of  cold  atoms  in  an  optical  lattice, 
which  is  periodic  potential  for  the  atoms  formed  by  interfering  several  laser  beams. 
Such  optical  lattices  allow  one  to  realize  analogous  physics  to  strongly  correlated 
electron  systems,  but  in  a  controlled  environment  with  much  less  noise. 

The  held  of  cold  gases  in  optical  lattices  is  by  now  a  well  developed  with  sev¬ 
eral  seminal  discoveries  including  the  observation  of  the  superfluid  to  Mott  insulator 
transition  in  the  Bose-Hubbard  model,  the  crossover  from  a  Bose-Einstcin  condensate 
(BEC)  to  a  Bardeen-Cooper-Schriefer  (BCS)  superfluid  and  the  quantum  phase  tran¬ 
sition  of  an  antiferromagnet  in  the  Ising  model  (Bloch  et  ah,  2012;  Grimm  et  ah,  2000; 
Simon  et  ah,  2011).  One  of  the  main  goals  of  quantum  simulation  in  optical  lattices 
is  to  realize  the  two  dimensional  Fermi-Hubbard  model  at  very  low  temperatures  to 
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determine  if  there  is  a  d-wave  superfluid  phase.  This  is  an  outstanding  question  in 
condensed  matter  physics  with  many  implications  for  high-Tc  superconductors  and 
other  strongly  correlated  systems.  To  explore  these  issues  improvements  to  existing 
optical  lattice  systems  must  be  made.  In  this  thesis  we  explore  a  novel  optical  lattice 
system  where  the  trapping  held  is  formed  from  the  scattered  light  off  an  array  of 
plasmonic  nanoparticles,  which  allows  one  to  increase  the  energy  scales  of  the  sys¬ 
tem  and  achieve  novel  long-range  interactions  via  the  electromagnetic  modes  of  the 
nanoparticles. 

1.2.3  Quantum  Nonlinear  Optical  Systems 

Recent  years  have  seen  many  breakthroughs  in  our  ability  manipulate  and  control 
light.  On  the  one  hand,  advances  in  materials  science  and  nanoscience  have  allowed 
the  design  of  devices  with  structure  well  below  the  wavelength  of  light.  Notable 
examples  include  photonic  crystals  in  dielectric  media  (Joannopoulos  et  al.,  2008), 
plasmonic  structures  in  metallic  systems  (Barnes  et  ah,  2008),  metamaterial  systems 
(Shalaev,  2007),  and  optomechanical  systems  (Marquardt  and  Girvin,  2009).  On  the 
other  hand,  the  growing  held  of  quantum  information  science  has  provided  a  new 
standard  for  controlling  the  quantum  properties  of  light,  as  well  as  a  host  of  novel 
platforms  to  achieve  this  control  (O’Brien  et  ah,  2009).  A  fundamental  challenge  for 
applications  of  these  systems  is  achieving  strong  interactions  between  photons.  The 
most  stringent  example  is  in  quantum  information,  where  one  requires  significant 
nonlinearities  at  the  level  of  a  single  quanta.  In  addition,  photonic  and  quantum 
optical  systems  offer  a  new  paradigm  in  theoretical  physics  in  that  they  are  funda- 
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mentally  non-equilibrium  systems.  Such  systems  provide  an  exciting  opportunity  to 
develop  new  technologies  in  both  the  classical  and  quantum  domains,  as  well  as  probe 
fundamental  questions  regarding  non-equilibrium  many-body  physics 

In  this  thesis  we  analyze  several  systems  where  it  is  possible  to  realize  nonlinear 
optical  effects  at  the  level  of  a  few  photons.  First  we  show  that  the  intrinsic  non¬ 
linearity  for  plasmons  in  graphene  nanostructures  is  strong  enough  that  the  material 
can  become  nonlinear  at  the  level  of  a  single  plasmon  as  illustrated  in  Figure  1  .  Such 
effects  occur  clue  the  subwavelength  confinement  of  the  plasmons  compared  to  free 
space,  which  significantly  enhances  the  electric  held  intensity  per  photon.  We  then 
go  on  to  look  at  nonlinear  effects  in  atomic  ensembles  where  the  long  coherence  times 
allow  one  to  store  the  photons  as  matter  for  long  times  to  achieve  large  interactions. 
We  consider  two  approaches,  one  based  on  cavity  quantum  electrodynamics  (QED) 
to  achieve  the  interactions  and  the  other  based  on  excitation  to  strongly  interacting 
Rydberg  states. 

1.3  Structure  of  Thesis 

Chapter  2  of  this  thesis  is  focused  on  double  dot  electron  spin  qubits.  We  show 
how  to  control  and  prepare  the  nuclear  spin  environment  of  the  electron  spins  through 
dynamic  nuclear  polarization.  In  Chapter  3  we  propose  and  analyze  a  novel  approach 
to  the  realization  of  high-density  optical  lattices  using  the  optical  potential  formed 
from  the  near  held  scattering  of  light  by  an  array  of  plasmonic  nanoparticles.  In 
chapters  4-6  we  consider  quantum  nonlinear  optical  systems.  In  Ch.  4  we  consider 
the  enhanced  nonlinearity  for  plasmons  in  graphene  nanostructures.  In  Ch.  5  we 
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theoretically  analyze  an  all-optical  single  photon  switch  using  slow  light  and  cavity 
QED.  Finally,  in  Ch.  6  we  consider  the  few  body  physics  of  strongly  interacting 
photons  in  Rydberg  systems  where  we  study  the  dynamics  of  two  and  three  body 
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Chapter  2 

Dynamic  Nuclear  Polarization  in 
Double  Quantum  Dots 

2.1  Introduction 

The  study  of  non-equilibrium  dynamics  of  nuclei  in  solids  has  a  long  history 
(Abragam  and  Goldman,  1978)  and  has  become  particularly  relevant  as  nanoscale 
engineering  and  improvements  in  control  allow  to  probe  mesoscopic  collections  of  nu¬ 
clear  spins  (Ynsa  et  ah,  2005;  Dixon  et  al.,  1997;  Salis  et  al.,  2001;  Ono  and  Tarucha, 
2004;  Koppens  et  al.,  2008;  Bracker  et  al.,  2005;  Lai  et  al.,  2006).  This  control  has 
direct  applicability  to  quantum  information  science,  where  nuclear  spins  are  often  a 
main  source  of  dephasing  (Hanson  et  al.,  2007).  The  goal  of  developing  an  under¬ 
standing  of  electronic  control  of  nuclei  is  to  circumvent  this  nuclear  dephasing  and  to 
turn  nuclear  spins  into  a  useful  resource  (Klauser  et  al.,  2008),  as  indicated  in  recent 
experiments  (Reilly  et  al.,  2008b;  Foletti  et  al.,  (2008,  2009;  Bluhm  et  al.,  2010,  2011; 
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Shulman  et  al.,  2012;  Frolov  et  al.,  2012). 

Double  quantum  dots  in  III-V  semicondnctors  can  be  operated  with  two  electrons 
coupled  to  approximately  104  to  106  nuclei  by  the  contact  hyperfine  interaction.  Re¬ 
peated  cycles  transitioning  from  the  electronic  singlet  to  triplet  states  can  be  used  to 
polarize  the  nuclear  spins;  electron  spin  flips  between  the  singlet  and  triplet  spaces 
occur  dne  to  the  difference  D  in  the  Overhauser  fields  on  the  two  dots  (Petta  et  ah, 
2008).  Early  experimental  (Reilly  et  ah,  2008b)  and  theoretical  (Ramon  and  Hu, 
2007;  Ribeiro  and  Burkard,  2009;  Yao  and  Luo,  2010;  Stopa  et  ah,  2010)  work  sug¬ 
gested  that  the  polarization  process  naturally  drove  the  projection  of  the  difference 
held  onto  the  magnetic  held  axis  Dz  to  zero.  However,  later  experiments  and  theory 
both  showed  that  the  polarization  is  naturally  accompanied  by  a  growth  in  Dz  and 
that  the  data  in  the  original  experiments  showing  a  suppression  in  Dz  was  likely  mis¬ 
interpreted  (Foletti  et  ah,  2009;  Gullans  et  ah,  2010).  Instead  the  results  are  more 
consistent  with  the  growth  of  a  large  Dz  accompanied  by  a  reduction  in  measurement 
contrast  between  singlet  and  triplet  states,  which  makes  it  appear  as  if  Dz  is  small 
(Barthel  et  ah,  2012). 

In  this  chapter  we  develop  a  model  to  describe  the  long  time  dynamics  of  the 
nuclear  spins  undergoing  adiabatic  pumping.  These  results  are  in  good  agreement 
with  the  experiments  described  above  (Foletti  et  ah,  2009;  Shulman  et  ah,  2012). 
The  main  conclusion  from  this  work  was  that  when  the  dots  are  different  sizes  the 
Overhauser  held  becomes  larger  in  the  smaller  dot;  thereby  resulting  in  large  differ¬ 
ence  fields.  In  the  present  work,  we  present  a  detailed  theoretical  analysis  of  these 
problems.  We  describe  the  theoretical  methods  developed  to  study  this  system,  in- 
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eluding  a  novel  method  for  efficient  simulation  of  semiclassical  central  spin  problems, 
and  detail  the  experimentally  relevant  polarization  phenomena  we  find  in  our  model. 
The  main  results  of  the  present  work  are  that  when  nuclear  spin  noise  is  included,  the 
more  detailed  theory  presented  here  agrees  with  the  results  of  Gullans  et  al.  (2010); 
however,  in  the  absence  of  nuclear  spin  noise,  states  with  D~  =  0  can  also  be  achieved 
for  certain  parameters. 

Our  theoretical  methods  are  based  on  a  semiclassical  description  of  the  nuclear 
spin  dynamics  in  which  the  nuclear  spins  are  grouped  into  small  sets,  each  homo¬ 
geneously  coupled  to  the  electron  spin  (Christ  et  ah,  2007).  The  nuclei  in  each  set 
may  be  treated  as  a  single  collective  spin  and  a  semiclassical  treatment  is  justified 
provided  the  number  of  spins  in  each  set  remains  large.  Increasing  the  number  of 
such  sets  improves  the  approximation  to  the  true  hyperhne  coupling.  More  formally, 
we  construct  a  systematic  approximation  to  the  true  hyperhne  coupling  in  terms  of  a 
reduced  set  of  M  coupling  constants.  For  the  optimal  choice  of  coupling  constants,  we 
rigorously  prove  that  our  approximation  reproduces  the  exact  semiclassical  time  dy¬ 
namics  to  within  a  fixed  error  for  a  time  that  increases  linearly  with  M.  For  large  M, 
this  allows  examination  of  the  long  timescales  relevant  for  polarization  experiments. 
This  approach  extends  previous  work  that  assumes  that  all  nuclei  on  a  given  dot 
have  equal  coupling  to  the  electron  spin  (Ramon  and  Hu,  2007;  Ribeiro  and  Burkard, 
2009;  Yao  and  Luo,  2010;  Stopa  et  ah,  2010;  Brataas  and  Rashba,  2011;  Rudner  and 
Levitov,  2012);  an  approach  which  often  incorrectly  predicts  rapid  saturation  of  the 
polarization.  Other  extensions  to  this  homogenous  coupling  model,  including  semi¬ 
classical  solutions  for  the  central  spin  (Brataas  and  Rashba,  2012;  Chen  et  al.,  2007; 
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Al-Hassanieh  et  al.,  2006;  Tsyplyatyev  and  Loss,  2011),  and  cluster  and  diagramatic 
expansion  techniques  for  short  time  non-equilibrium  behavior  (Witzcl  and  Sarma, 
2008;  Wang  et  ah,  2006;  Coish  and  Loss,  2004)  do  not  explore  the  wide  range  of  time 
scales  or  relevant  physics  for  the  double  dot  case. 

Our  results  can  be  broken  up  into  two  distinct  cases  depending  on  whether  or  not 
the  dots  are  identical.  When  the  dots  are  different  sizes,  then  the  hyperhne  coupling, 
which  scales  inversely  with  the  volume,  is  larger  on  the  smaller  dot  and  we  find 
that  the  Overhauser  field  grows  preferentially  on  the  smaller  dot  as  the  polarization 
increases.  This  preferential  growth  results  in  a  large  Overhauser  difference  field  Dz. 
For  two  dots  with  a  difference  in  volume  of  less  than  ~  20%  we  find  a  rich  and  complex 
phase  diagram  for  the  nuclear  spin  dynamics,  which  can  be  broken  into  two  distinct 
regimes.  The  first  regime  occurs  with  large  external  magnetic  fields  or  short  cycle 
times.  In  this  regime  the  system  saturates  without  significant  polarization  because  the 
perpendicular  components  of  D  rapidly  approach  zero  and  spin  flips  are  suppressed; 
the  system  approaches  a  semiclassical  “dark  state.”  This  occurs  with  no  statistical 
change  in  the  distribution  of  Dz.  The  second  regime  occurs  in  the  limit  of  smaller 
magnetic  fields  or  slower  cycle  times.  In  this  regime,  the  dynamics  are  sensitive  to 
the  inclusion  of  nuclear  spin  noise.  In  the  absence  of  nuclear  spin  noise  we  find  one 
potential  end  state  of  polarization  is  a  “zero  state”  in  which  all  components  of  D  — Y  0. 
In  this  state  the  singlet  and  triplet  electronic  subspaces  are  completely  decoupled  and 
spin  flips  no  longer  occur.  Simultaneously,  though,  there  are  instabilities  leading  to 
the  growth  of  large  Overhauser  difference  fields.  Crucially,  when  even  a  small  amount 
of  nuclear  spin  noise  is  added  the  zero  states  strongly  destabilize  and  the  system 
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generically  becomes  unstable  to  the  growth  of  large  difference  fields  as  shown  by 
Gullans  et  ah  (2010). 

These  results  provide  a  clear  picture  of  the  polarization  dynamics  in  such  double 
quantum  dot  systems  and  will  be  a  useful  guide  to  future  experiments  aimed  at  more 
precise  control  of  the  nuclear  spins.  Although  the  paper  is  specific  to  double  quantum 
dots  in  GaAs,  many  of  the  results  and  theoretical  methods  extend  to  other  central 
spin  systems  under  investigation  (Takahashi  et  al.,  2011;  Sun  et  al.,  2012;  Hogele 
et  ah,  2012).  More  generally,  this  work  is  of  fundamental  interest  as  we  explore 
the  dynamics  of  an  interacting,  many-body  system  when  it  is  far  from  equilibrium 
(Urbaszek  et  ah,  2013). 

The  paper  is  organized  as  follows.  In  section  2  we  define  the  Hamiltonian  for  the 
double  dot  system  and  introduce  the  polarization  cycle.  In  section  3  we  systematically 
derive  a  semiclassical  model  for  the  nuclear  spins  starting  from  the  coarse-grained 
evolution  of  the  nuclear  spin  density  matrix.  In  section  4  we  present  our  results  for 
identical  and  unequal  dots  in  the  presence  and  absence  of  nuclear  spin  noise.  In 
appendix  A.l  we  provide  a  summary  of  the  parameters  used  in  our  simulations.  In 
appendix  A. 2  we  describe  our  approach  to  coarse  graining  the  electron  wave  function 
and  provide  rigorous  bounds  on  the  error  in  time  evolution  due  to  the  course  graining. 
In  appendix  A. 3  we  extend  our  simulations  to  the  case  of  multiple  nuclear  species  and 
find  qualitatively  the  same  results  as  for  a  single  species. 
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Figure  2.1:  a)  The  Overhauser  field  in  each  dot  gives  rise  to  sum  and  difference  fields 
which  are  relevant  for  the  double  dot  system,  b)  Schematic  of  two-electron  energy 
levels  as  a  function  of  detuning  e  between  (1,1)  and  (0,2)  charge  states.  Arrows 
indicate  adiabatic  sweep  through  avoided  crossing  (pink)  and  rapid  sweep  back  to 
(0,2)  with  reload  (green),  c)  Spin-flip  pathways  between  the  s  and  T+  states  as  the 
exchange  energy  J(e)  is  swept  through  the  crossing,  showing  the  nuclear  operators 
involved  in  each  path.  Each  pathway  is  a  term  in  D_  in  Eq.  2.2. 


2.2  Setup 

For  a  double  quantum  dot  with  two  electrons,  we  can  write  the  Hamiltonian  for  the 
lowest  energy  (1,1)  and  (0,2)  electron  states,  where  (n,m.)  indicates  n  ( m )  electrons 
in  the  left  (right)  dot.  To  model  nuclear  polarization,  we  first  derive  an  effective  two- 
level  Hamiltonian  to  describe  the  system  near  the  crossing  of  the  singlet  s  and  lowest 
energy  triplet  state,  T+,  of  this  two-electron  system,  then  solve  the  time  dynamics. 
Dynamic  nuclear  polarization  (DNP)  experiments  operate  near  this  crossing,  typically 
with  an  adiabatic  sweep  of  the  difference  in  the  dots  electric  potential  through  the 
s-T+  degeneracy  (Fig.  2.1a),  followed  by  a  non-adiabatic  return  to  (0,2)  and  reset  of 
the  electronic  state  via  coupling  to  leads. 

If  ipd(r)  is  the  single-particle  envelope  wave  function  on  dot  d  =  l, r  (for  the 
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left,  right  dot),  the  effective  hyperhne  coupling  for  the  nuclear  spin  at  rkd  is  gkd  = 
ahfVo\ipd(rkd)\2  where  ahf  is  the  hyperhne  coupling  constant,  and  no  is  the  volume  per 
nuclear  spin.  We  introduce  two  collective  nuclear  spin  operators  to  denote  the  Over- 
hauser  fields  in  the  left  (L)  and  right  (R)  dots,  L  =  gkdki  and  R  =  gfcrIfcr,  and 
further  dehne  S  —  (L  +  R)/ 2,  D  =  (L  —  R)/ 2,  where  I kd  is  the  angular  momentum 
of  the  kth  nucleus  on  dot  d.  The  rms  Overhauser  energy  in  the  infinite  temperature 
ensemble  is  Dd  =  +  l)/3)1^2  where  /  is  the  magnitude  of  each  nuclear  spin. 

We  dehne  fl  =  \J (fi2  +  Dj)/2,  and  work  in  energy  and  magnetic  held  units  such  that 
Q  =  —  2-^  =  1,  where  g*  is  the  electron  ehective  g-factor  and  /jlb  is  the  Bohr  magne¬ 
ton.  In  the  basis  (|s) ,  | T+) ,  |T0) ,  |T_)},  where  the  Tm  are  the  (1, 1)  triplet  states  and 
s  is  the  (1,  l)-(0,2)  hybridized  singlet  state,  the  Hamiltonian  is  (Taylor  et  al.,  2007) 


-J(e) 

vD+ 

—  \/2vDz 

—vD- 

v  D _ 

—  Bex  t  +  Sz 

S-/V 2 

0 

—  y/2  vDz 

S+/V2 

0 

S-/V 2 

—vD+ 

0 

S+/V 2 

Text  -  5; 

where  D±  =  Dx  ±  iDy  and  similarly  for  S±,  Bext  is  an  external  magnetic  held, 
v  =  v{e )  =  cos6)(e)/v/2,  and  cos 9{e)  is  the  overlap  of  the  (1,1)  singlet  state  with 
the  (l,l)-(0,2)  hybridized  singlet  state  s).  The  parameters  cos 0(e)  and  J(s),  the 
splitting  between  s  and  To,  are  both  functions  of  the  energy  difference  £  between 
the  (1,1)  and  (0,2)  charge  states.  Here  the  nuclear  spin  variables  refer  to  the  full 
quantum  mechanical  operators  on  the  nuclear  spin  space.  In  appendix  A. 3  we  will 
consider  the  case  of  multiple  nuclear  species,  but  for  now  we  consider  the  nuclei  to 
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be  spin-3/2  of  a  single  species,  in  a  frame  rotating  at  the  nuclear  Larmor  frequency. 

Assuming  that  J,  Bext  A>  fl,  we  perform  a  formal  expansion  in  the  inverse  electron 
Zeeman  energy  operator  m  =  D/(Bext  —  Sz  +  irj)  where  rj  >  0  is  infinitesimal.  We 
apply  a  unitary  transformation  that  rotates  the  quantization  axis  of  the  triplet  states 
to  align  with  Bext  —  S  and  find  the  Hamiltonian  for  the  (|s) ,  |X+)}  subspace  to  first 
order  in  J-1,  m: 


HcS  = 


—  J{s)  +  hs 
K  v(e)D_ 


v(e)D+ 
—Bext  +  hr j 


(2.1) 


where  the  effect  of  coupling  to  the  higher  energy  states  |T0)  and  \T_)  enters  as 


hs  =  -—D\DZ  -  D_ - — - - D+ , 

J  Z  J  +  Bext-Sz  +’ 

hT  =  SZ-  i{S-S+m  +  mS.S+), 

D_  =  D-  +  rhS-bz  —  fhS2_mb+  —  |mS_5+mD_t 
b  z  =  D  z  —  -(S+mb_  +  S_rhb+). 


Of  particular  interest  is  that  the  off-diagonal  term,  which  produces  nuclear  polariza¬ 
tion,  vanishes  in  the  semiclassical  limit  of  (D)  — y  0,  i.e. ,  in  the  zero  states. 


2.3  Model 

We  develop  a  model  for  the  evolution  of  the  nuclear  spin  density  matrix  after  one 
pair  of  electrons  has  cycled  through  the  system.  We  approximate  the  sweep  through 
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the  |s)-|T+)  degeneracy  as  a  Landau-Zener  process,  which  we  solve  approximately 
for  the  effect  on  the  nuclear  system.  By  coarse-graining  this  evolution  over  a  cycle 
we  derive  a  master  equation  for  the  nuclear  spins.  Finally,  we  add  the  effects  of 
nuclear  dipole-dipole  interactions  and  qnadrupole  splittings  phenomenologically.  The 
derivation  presented  here  is  complementary  to  that  of  Gullans  et  ah  (2010)  and  results 
in  the  same  equations  of  motion. 

The  electron  system  is  prepared  in  |s)  at  large  negative  t  =  —  T/2,  where  T  is  the 
total  cycle  time.  We  identify  the  (nuclear  spin)  eigenstates  of  the  operator  D+D_, 
labeled  \Dj_)  with  eigenvalues  D\.  Since  the  components  of  hs  and  hr  that  do  not 
commute  with  D+D_  are  perturbatively  small  in  m0  and  1/J,  we  approximate  them 
by  keeping  only  the  diagonal  components  in  the  two-level-system  subspace,  sending 
hs  — *  (D_ l|  hs  \D±)  and  hT  — »  (D'±  \D'±)  where  | D'±)  =  D~^D_  \D±).  In  this  limit, 

the  off-diagonal  part  of  Hc g  in  Eq.  2.1  produces  standard  Landau-Zener  behavior, 
while  the  diagonal  components  of  Heg  are  simply  phases  picked  np  by  the  nuclei, 
depending  on  which  electronic  state  is  occupied.  For  initial  state  |T0)  =  |s)  <8>  \D±), 
the  crossing  either  leaves  the  electronic  state  unchanged  or  flips  an  electron  and 
nuclear  spin  to  the  state  |T+)  (8)  | D'±).  We  note  that  \D'±)  is  an  eigenstate  of  D~D+ 
with  eigenvalue  D'j_.  The  problem  is  now  reduced  to  finding  Landau-Zener  solutions 
for  each  independent  two-level  system  (|s)®  \D±),  \T+)  <g)  |Z)^_) }.  We  model  the  actual 
sweep  of  £  by  a  linear  sweep  of  J  so  J(t)  =  —2 /32t+Bext,  where  (3  =  \J |  \dJ{e)/dt\  |t=o- 
We  take  v(e)  to  be  constant,  valid  in  the  limit  of  large  tunnel  coupling,  and  assume 
/ 3  <C  Bext  to  ensure  the  applicability  of  Eq.  2.1.  For  moderate  magnetic  fields  v(e)  ~ 
l/\/2,  but  it  decreases  at  large  magnetic  fields  as  the  (l,l)-(0,2)  hybridized  singlet 
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state  has  a  smaller  overlap  with  (1,1)  at  the  s-T+  crossing. 

After  one  cycle,  |^0)  evolves  into  |^i)  =  cs  |s)  ®  \D±)  +  ct\T+)  ®  \ D'±).  For 
f32T  >>  1,  the  standard  Landau-Zener  formula  gives  the  flip  probability  as  pf  = 
1  —  exp(— 27 tu2),  where  oj  =  v{D±)/j3,  and 


cs  =  y/l-Pfexp(-i^>s),  cT  =  y/pj  exp(—i(/>T) 

fT/2 

4>s  ~  /  hsdt 


(2.3) 


4>t 


’  —T/2 
ft  o 

’  —T/2 


hsdt  +  {T/2  —  to)  hr  +  4>ad{^), 


where  the  crossing  occurs  at  a  time  to  ~  Sz/ f32.  We  include  in  the  phase  picked 
up  by  following  the  adiabat,  4>ad-  We  approximate  (Pad  by  interpolating  between  the 
limits  u  =  v(Dj_)//3  — »  0  and  uj  — >  oo,  giving  (Vitanov  and  Garraway,  1996) 


4>ad 


27TC02  +  Pf 


1  —  2n  +  log 


5 


where  r  =  Tfd/2.  More  accurate  approximations  can  easily  be  taken  into  account 
within  our  formalism;  however  we  find  such  corrections  have  a  negligible  effect  on  the 
long  term  polarization  dynamics  because  the  polarization  process  rapidly  drives  u  to 
small  values. 

We  move  from  the  independent  two-level  systems  to  the  general  case  by  noting 
that  the  components  of  |\P)  depend  only  on  the  eigenvalue  Dj_  and  on  the  polarization 
Sz  (which  we  approximate  as  commuting).  Since  the  eigenstates  of  D+D _  form 
a  complete  basis  for  the  nuclear  spin  states  we  can  define  the  complete  operator 
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pf  =  J2D±Pf(D±)  \Dj_)  (Dj_\,  and  similarly  for  05,  0T.  The  nuclear  spin  density 
matrix  after  each  cycle  is  given  by  tracing  over  the  electronic  states.  The  nuclear 
density  matrix  evolution  is  then 


where  pn  is  the  nuclear  density  matrix  after  n  cycles. 

Rather  than  solve  for  the  exact  dynamics  of  the  nuclear  density  matrix-still  an 
intractably  hard  computational  problem  for  any  reasonable  number  of  nuclear  spins- 
we  instead  adopt  an  approximate  solution  to  the  problem  using  the  P-representation 
for  the  density  matrix  as  an  integral  over  products  of  spin  coherent  states.  From  the 
thermal  distribution,  we  choose  such  a  spin  coherent  state  and  evolve  it,  where  we 
interpret  expectation  values  (...)  as  being  taken  in  that  state.  The  ensemble  of  such 
trajectories  represents  the  physical  system  (Al-Hassanieh  et  al.,  2006). 

We  organize  this  calculation  by  noting  that  the  components  of  the  Landau- Zener 
model  (05,  (j)T,pf,  D±)  are  only  functions  of  L  and  R.  A  spin  coherent  state  is  entirely 
described  by  its  expectation  values  i u  =  (1^).  For  the  kth  spin  on  the  left  dot,  we 
expand  the  discrete  time  difference  {lki)n  —  (Ifc0n_1  after  n  and  n  —  1  cycles  in  the 
small  parameter  ,  giving  an  evolution  equation 

3 

=  9ki  y:  Pi,n  (^[9gklLfj,,lki]j  =  gkiR  x  hi,  (2-4) 

/^=i 


19 


38 


Chapter  2:  Dynamic  Nuclear  Polarization  in  Double  Quantum  Dots 


with 


Pl  =  7p 


(l-pf){Vi4>s)  +  (P/)(V;0r)  -  Im(7,) 


where  V/  =  {dLx,dLy,dLz)  and 

=  (2’5) 

and  similarly  for  ifcr,  P,.,  and  qy,  with  L  replaced  by  R.  The  factorization  of  expec¬ 
tation  values  is  a  natural  consequence  of  our  spin-coherent  state  approximation,  as 
it  explicitly  prevents  entanglement  between  spins.  Thus  we  have  an  effective,  semi- 
classical  picture  of  nuclear  spins  precessing  and  being  polarized  by  their  interaction 
with  the  electron  spin,  integrated  over  one  cycle. 

We  approximate  the  electron  wavefunction  as  a  piecewise-flat  function  with  M 
levels,  which  we  refer  to  as  the  annular  approximation,  as  illustrated  in  Fig.  2.2a. 
Each  annulus  defines  I nd  =  Yhk&n  hw,  where  the  sum  is  over  all  nuclei  with  the  same 
hyperhne  coupling  to  the  electron.  Since  r/fc  is  identical  for  all  k  G  n,  we  can  simply 
replace  i kd,  with  lnd  in  Eq.  2.4.  Furthermore,  1%  is  a  conserved  quantity,  so  we  can 
study  the  evolution  of  M  <C  N  spins  in  a  reduced  Hilbert  space.  The  typical  size  of 
In  is  ~  y/W/M  3>  1,  which  allows  us  to  replace  the  spin-coherent  states  used  above 
with  semi-classical  spins,  and  makes  taking  expectation  values  straightforward:  all 
quantum  operators  can  be  replaced  by  their  expectation  values  directly.  The  annular 
approximation  should  correctly  describe  the  nuclear  dynamics  for  a  time  scale  given 
by  the  inverse  of  the  difference  between  the  r/fc  of  adjacent  annuli. 
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To  illustrate,  to  first  order  in  m0  =  £>ext,  for  d  =  l,r, 


Pd  =pf\ (A +z  -  A0S±)  +  mor0|^|^  x  D 

r  fi2 

+  TRpfVd(j)AD  =F  [r0^Im  (7 1  ~  lr) 

+  (1  —  pf\/2)(AoDzz  +  A_Dj_) 


(2.6) 


where  the  top  sign  applies  for  d  —  /,  =  (Tb,  Dy,  0),  =  (Sx,  Sy,  0),  A  =  1  —  2 t0/T 

gives  the  shift  in  the  location  of  the  crossing,  and  A0,  A_,  A+,  A0,  Tr,  and  T0  are 
constants  depending  on  the  details  of  the  pulse  cycle  (see  below).  We  have  replaced 
operators  by  their  expectation  values  and  removed  the  angle  brackets  since  we  are  now 
in  the  semiclassical  limit.  To  leading  order  in  mo,  Im  (7;  —  7r)  =  2(D  x  z)pf/D\.  It  is 
clear  from  Eqs.  2. 4-2. 6  that  all  dynamics  stop  in  the  zero  states  with  D  =  0,  consistent 
with  the  idea  that  true  saturation  of  polarization  requires  that  all  components  of  D 
be  small.  We  will  focus  on  the  stability  of  such  states  in  various  parameter  regimes. 
The  equations  of  motion  in  Gullans  et  al.  (2010)  are  found  from  Eq.  2.6  by  including 
only  the  lowest  order  in  DT  and  0//3,  which  is  the  limit  of  fast  cycles  and  small  spin 
flip  probability  per  cycle,  respectively. 

First  we  outline  the  meanings  of  the  parameters  in  the  model.  As  indicated 
schematically  in  Fig.  2.2b,  the  To  term  originates  in  the  hyperhne  flip-flop,  the  Ao 
and  A_  terms  are  the  off-resonant  effects  of  coupling  from  the  singlet  state  to  the  To 
and  T_  states,  respectively,  A0  comes  from  coupling  between  the  T+  and  T0  states, 
and  A+  comes  from  Knight  shifts  due  to  occupation  of  the  T+  state.  To  leading 
order  in  mo,  for  a  pulse  sequence  consisting  of  only  the  Landau-Zener  sweep,  with 
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Figure  2.2:  a)  Independent  Random  Variable  Annular  Approximation  (IRVAA)  to 
the  electron  wavefunction  in  the  double  dot.  b)  Key  processes  contributing  to  Eq.  2.6. 


instantaneous  eject  and  reload,  the  parameters  have  values 


A+  =  1/4, 

7-1  2? rv2fc 

r°  ~ 


m0, 


A0  —  m0 /A 

r  R  =  fe 


where  fc  —  1/T  is  the  cycle  frequency  and  (.)c  indicates  an  average  taken  over  a  full 
cycle;  these  values  can  be  modified  readily  by  changing  the  details  of  the  pulse  cycle, 
while  leaving  the  Landau-Zener  portion  unchanged.  In  Appendix  A.l  we  provide  a 
reference  for  all  parameters  used  in  the  simulations. 

Equation  2.4  is  a  good  approximation  of  the  nuclear  dynamics  over  a  few  DNP 
cycles  because  other  nuclear  processes  are  slow  compared  to  a  typical  experimental 
cycle  (~10-100  ns  (Reilly  et  al.,  2008b)).  However,  the  full  DNP  may  last  millions  of 
cycles  at  which  point  these  other  nuclear  processes  become  important.  Apart  from 
Larmor  precession,  which  is  only  relevant  for  the  case  of  multiple  nuclear  species 
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considered  in  Appendix  A. 3,  nuclear  qnadrnpole  splittings  and  nuclear  dipole-dipole 
interactions  are  the  dominant  processes.  They  become  relevant  on  a  timescale  of  a 
few  hundred  microseconds  in  these  systems  (Taylor  et  al.,  2007).  We  include  them 
in  our  model  phenomenologically  by  adding  a  fluctuating  magnetic  field  hkd(t )  in  the 
^-direction  at  each  site  (the  transverse  terms  are  strongly  suppressed  by  the  external 
held),  such  that 

“77  9kd-Pd  X  ikd  bn  ^ kd  %  X  ikd  (2-7) 

at 

where  yn  is  the  nuclear  gyromagnetic  ratio.  We  further  assume  that  the  this  held 
can  be  treated  as  noise  and  characterized  by  a  Gaussian,  uncorrelated  white  noise 
spectrum 

lliKdi^K'd'i^n  =  2 r}S(t  ~  t')8kk'8ddl  (2.8) 

where  (-)n  are  averages  over  the  noise  (Reilly  et  ah,  2008a). 

2.4  Results 

The  polarization  dynamics  display  three  characteristic  behaviors:  growth  of  large 
difference  helds,  saturation  in  nuclear  dark  states  defined  by  D±  =  0,  and  preparation 
in  zero  states  D  —  0  which  are  global  hxed  points  of  the  nuclear  dynamics  in  the 
absence  of  noise.  In  Gullans  et  ah  (2010)  this  system  was  studied  in  a  restricted  model 
focusing  on  the  case  where  noise  was  present.  Therein  it  was  found  that  when  the  two 
dots  have  different  hyperhne  couplings  the  system  generically  grows  large  difference 
holds,  while  for  identical  dots,  depending  on  parameters,  the  system  is  either  unstable 
to  the  growth  of  large  difference  helds  or  saturates  in  dark  states;  however,  the  zero 
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states  were  not  found  to  be  a  relevant  steady  state  in  any  parameter  regime.  In 
the  present  work  we  focus  on  extending  the  results  of  Gullans  et  al.  (2010)  to  a 
larger,  more  experimentally  relevant,  parameter  regime  by  using  equations  of  motion 
correct  to  second  order  in  mo  with  a  more  complete  model  of  the  Landau-Zener 
sweep  as  described  in  the  previous  section.  In  addition,  we  consider  the  nuclear 
dynamics  in  the  absence  of  noise.  We  also  present  the  full  analytical  calculations 
which  were  omitted  from  Gullans  et  al.  (2010).  In  all  physical  parameter  regimes  we 
find  qualitatively  consistent  results  with  Gullans  et  al.  (2010);  however,  for  a  limited, 
unphysical  parameter  regime  we  do  find  solutions  to  the  equations  of  motion  in  the 
absence  of  noise  where  the  zero  state  is  uniformly  reached  starting  from  a  completely 
uncorrelated  nuclear  spin  ensemble. 

The  simulations  shown  below  were  performed  with  the  equations  of  motion  correct 
to  second  order  in  m0  with  ^(r)  a  2D  Gaussian.  Taking  v2  «  1/2,  we  estimate  that 
for  experiments  performed  with  Bext  =  10  mT  with  T  =  25  ns  (Reilly  et  ah,  2008b), 
mo  ~  0.18,  T0  ~  0.20,  but  the  A  and  A  terms  depend  on  the  rest  of  the  cycle.  In 
each  of  the  simulations,  we  choose  initial  magnitudes  and  directions  of  the  spins  In 
by  a  procedure  equivalent  to  choosing  initial  directions  for  each  of  the  Nn  spin-3/2 
nuclei  in  the  nth  annulus  and  evaluating  I„  =  i k  explicitly  (see  Appendix  A. 2). 
The  relationship  between  simulation  time  and  laboratory  time  depends  on  the  details 
of  the  pulse  cycle,  including  pauses  and  reloads  not  considered  explicitly  here,  but 
simulation  time  is  roughly  in  units  of  where  gmax  ~  2f I2/a/,/  is  the  largest  value 
of  gfc,  so  t  =  400  is  approximately  10  ms. 

To  organize  our  results  we  recall  the  phase  diagram  for  identical  dots  and  the 
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0  0.2  0.4  0.6  0.8 


A_/A0 

Figure  2.3:  Phase  diagram  for  the  simplified  model  presented  in  Gullans  et  al.  (2010). 
At  each  value  of  parameters,  twenty  runs  were  started  with  Dz  =  —2,  Sz  =  —10,  and 
all  other  components  chosen  randomly  according  to  the  infinite  temperature  ensemble. 
The  colorscale  indicates  how  many  of  those  runs  ended  with  \DZ\  increased.  The  dark 
region  is  of  saturation  and  the  light  region  is  of  instability.  The  dashed  line  shows 
the  prediction  of  the  simple  model  of  Eq.  2.30,  which  captures  the  phase  boundary, 
especially  at  low  A_/Aq.  For  parameters  used,  see  Table  I. 
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Figure  2.4:  Phase  diagram  as  in  Fig.  2.3,  except  with  varying  external  magnetic 
field  and  without  any  noise  added.  The  parameters  were  scaled  with  mo  as  shown  in 
Section  III.  There  is  a  clear  boundary  between  saturation  at  large  T0  and  instability 
at  lower  values  of  r0,  with  appropriately  large  values  of  A0  and  A_.  See  Table  I  for 
parameters.  The  symbols  ’x’  and  ’o’  mark  the  parameters  used  for  Fig.  2.7  below. 
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simplified  model  derived  in  Gullans  et  al.  (2010)  in  which  the  only  non-zero  parameters 
are  Ao,_,  To  and  rj,  which  corresponds  to  the  limit  of  large  magnetic  fields  and  fast 
sweeps  including  nuclear  spin  noise.  To  obtain  the  phase  diagram  we  consider  for 
each  set  of  parameters  whether  the  system  supports  self-consistent  growth  of  \DZ\ 
starting  from  large  values  of  \DZ\  and  \SZ\.  This  approach  avoids  complications  with 
the  metastability  of  zero  states  discussed  later.  Such  simulations  produce  the  phase 
diagram  in  Fig.  2b  of  Gullans  et  al.  (2010),  which  is  reproduced  in  Fig.  2.3  with  the 
full  data  presented.  From  this  figure  it  is  clear  that  we  can  separate  the  dynamics  into 
two  regimes  depending  on  parameters.  For  large  ratios  of  r0/A0,  which  corresponds 
to  large  magnetic  fields  or  strong  pumping  the  system  quickly  saturates  with  no 
growth  of  large  difference  fields.  For  small  ratios  there  is  an  instability  towards  large 
difference  fields.  In  the  first  section  we  explore  the  dynamics  in  the  absence  of  noise 
for  identical  dots  with  all  parameters  included.  In  the  second  section  we  include 
nuclear  spin  noise  and  asymmetry  in  the  dot  sizes. 

2.4.1  Noise  Free  Nuclear  Spins 

From  the  general  arguments  given  in  the  introduction  it  is  clear  that  when  the 
dots  have  different  hyperhne  couplings  the  system  naturally  grows  a  large  difference 
held.  Furthermore,  in  Gullans  et  al.  (2010)  it  was  shown  that  even  identical  dots 
display  similar  behavior  in  the  presence  of  noise.  In  this  section  we  analyze  the 
case  of  identical  dots  in  the  absence  of  noise  to  better  understand  the  role  of  the 
coherent  nuclear  dynamics.  We  begin  by  deriving  a  phase  diagram  analogous  to 
the  one  obtained  in  the  presence  of  noise  except  we  now  look  in  the  space  of  the 
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Figure  2.5:  a-b)  Simulations  corresponding  to  the  saturation  region  of  the  phase 
diagram.  The  solid  lines  are  the  median  values  of  |DZ|  (a),  Sz  (a)  and  D±  (b)  at  each 
time  step  in  an  ensemble  of  1000  trajectories.  In  all  plots  shaded  regions  show  the  84th 
and  16th  percentiles,  c)  Simulations  showing  growth  of  (\DZ\)  with  the  time  shifted 
for  each  trajectory  so  that  its  maximum  \DZ\  occurs  at  time  zero.  Bottom  shows  the 
median  value  of  S2  (Sz)e  at  each  time  step  in  an  ensemble  of  1000  trajectories.  In  the 
middle  is  similar  (|D,|)e.  Thin  red  line  is  a  single  trajectory.  The  curve  at  top  shows 
the  fraction  of  trajectories  contributing  to  the  ensemble  at  each  time;  this  increases 
with  time  because  some  trajectories  reach  their  maximum  D.  much  later  than  others 
while  the  simulation  time  is  fixed  for  each  trajectory.  4.5%  of  the  trajectories,  which 
do  not  show  this  peak  in  \DZ\,  are  not  included.  Approximately  10%  of  the  trajectories 
show  behavior  similar  to  that  shown  in  the  thin  red  line,  where  \DZ\  is  reduced  initially 
and  then  goes  unstable  to  large  \DZ\.  d)  Mean  of  the  maximum  value  of  \DZ\  reached 
on  each  trajectory  for  the  same  parameters  as  in  (c)  (open  circles)  except  M  varied 
between  20  and  160,  with  5000  trajectories  per  point.  Closed  circles  show  similar 
results  with  m0  =  0.05,  t  —  4  and  all  other  parameters  scaled  appropriately.  The 
physical  system  has  M  — >  N  «  106,  so  we  interpret  this  as  an  instability  to  large  \DZ\, 
which  is  supported  by  simulations  including  transverse  noise  (see  section  2.4.2).  e) 
With  different  parameters,  simulations  showing  reduction  of  (\DZ\)  plotted  as  in  (c) 
without  the  time  shift.  For  these  parameters,  the  trajectories  have  \DZ\  — y  0  quickly, 
without  time  for  strong  polarization. 
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experimentally  accessible  parameters  cycle  rate  fc  and  inverse  magnetic  field  mo-  The 
results  are  shown  in  Fig.  2.4  where  we  see  the  same  qualitative  behavior  as  shown  in 
Fig.  2.3.  However,  the  dynamics  are  much  richer  than  indicated  by  this  simple  phase 
diagram.  In  the  following  subsections  we  give  examples  of  what  happens  to  a  nuclear 
spin  ensemble  starting  from  equilibrium  for  different  parameters  and  regions  of  the 
phase  diagram. 

Before  proceeding,  however,  we  note  that  in  the  absence  of  noise  the  inhomogeneity 
of  the  electron  wavefunction  plays  a  crucial  role.  This  is  because  weak  inhomogeneity 
is  equivalent  to  choosing  the  number  of  annuli  M  to  be  small  and  in  this  case  the 
system  moves  rapidly  to  its  maximally  polarized  state,  with  In  «  —Inz  for  all  n. 
Dynamics  completely  cease  in  this  state,  as  can  clearly  be  seen  from  Eq.  2.4,  despite 
the  fact  that  this  state  does  not  correspond  to  all  of  the  nuclei  being  polarized,  which 
would  also  require  In  =  3Nn/2.  On  the  other  hand,  for  strong  inhomogeneity,  or 
large  M,  when  the  system  is  not  fully  polarized  other  terms  in  Pci  compete  with  the 
polarization  saturation  and  sustain  the  dynamics.  (Christ  et  al.,  2007) 

Polarization  Saturation 

When  the  magnetic  held  is  large  or  the  cycle  rate  is  fast  (i.e.,  A0  <C  r0),  the 
system  rapidly  moves  toward  dark  states  (i.e.,  states  with  D±  =  0),  sending  pf  — >  0 
without  statistical  change  in  the  distribution  of  Dz,  as  shown  in  Fig.  2.5a.  This  limit 
is  additionally  characterized  by  only  a  small  change  in  nuclear  polarization  as  seen 
in  Fig.  2.5b.  When  the  effects  of  the  |s)-|T0)  coupling  are  important  (i.e.,  A0  m  r0), 
the  A0  term  in  Eq.  2.6  causes  D±  to  increase,  “rebrightening”  the  D±  ze  0  dark 


29 


48 


Chapter  2:  Dynamic  Nuclear  Polarization  in  Double  Quantum  Dots 

states  and  allowing  dynamics  to  continue.  Coupling  from  the  singlet  to  the  T0  state 
is  an  essential  ingredient  in  all  of  the  effects  discussed  below.  When  Ao  is  significant, 
dynamics  only  stop  near  zero  states  with  D  =  0. 

Growth  of  Difference  Fields 

Second,  we  observe  the  growth  of  large  Overhauser  fields.  We  consider  a  proto¬ 
typical  pulse  sequence  motivated  by  experiments  with  moderate/large  magnetic  held, 
m0  =  0.01  In  this  case,  over  95%  of  the  trajectories  display  a  growth  in  \DZ\,  as 
shown  in  Fig.  2.5c.  We  observe  this  behavior  over  a  range  of  experimentally  accessi¬ 
ble  magnetic  fields  and  cycle  frequencies.  This  increase  in  \DZ\  indicates  that  the  spin 
hips  are  occurring  predominantly  in  one  dot.  We  interpret  these  results  as  showing 
a  continuing  increase  of  \DZ\,  where  the  peak  of  \Dz(t)\  is  an  artifact  of  the  annular 
approximation.  Near  the  peak,  many  of  the  annular  spins  artificially  reach  their  max¬ 
imal  polarization,  at  which  point  they  should  be  broken  into  more  annuli.  Similar 
trajectories  with  different  M  show  the  maximum  value  of  \DZ\  increasing  with  M 
(Fig.  2.5d).  The  physical  cause  of  this  increase  in  \DZ\  is  not  clear,  but  it  is  associ¬ 
ated  with  both  A0/T0  and  A+/T0  being  sufficiently  large.  When  nuclear  spin  noise 
is  included,  the  growth  in  (| Dz\)e  continues  (Gullans  et  ah,  2010).  This  could  be  the 
same  phenomenon  as  seen  by  Foletti  et  al.  (2009),  with  transverse  dephasing  helping 
to  produce  the  large  \DZ\  ~  Bext  of  that  work,  though  unequal  dot  sizes  could  also 
produce  that  effect  (Gullans  et  ah,  2010). 
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Zero  States 

For  moderate  to  small  magnetic  fields,  when  Ao  ~  To,  two  different  characteristic 
behaviors  of  particular  note  are  observed.  First,  in  the  physical  parameter  regimes, 
which  do  not  display  general  motion  to  zero  states,  the  zero  states  are  still  important 
for  the  dynamics  as  they  are  a  metastable  state.  That  is,  many  trajectories  spend  a 
long  time  with  \DZ\  near  zero  before  escaping  away  to  large  \DZ\.  This  phenomenon 
is  shown  in  the  individual  trajectory  (thin  red  line)  of  Fig.  2.5c. 

Second,  for  parameters  in  our  model  which  are  not  experimentally  accessible  there 
is  a  mechanism  that  gives  rise  to  attraction  towards  zero  states.  This  is  illustrated 
in  Fig.  2.5e,  where  we  show  an  ensemble  of  trajectories  in  which  D  rapidly  reduces 
toward  zero.  For  the  parameters  of  Fig.  2.5e,  the  standard  deviation  of  Dz  was 
reduced  by  a  factor  of  28.  We  remark  that  as  D  — >  0,  the  singlet  state  ceases  mixing 
with  the  triplets  and  nuclear  spin  dynamics  stop.  Until  something  (outside  this  model, 
such  as  nuclear  dipole-dipole  coupling)  restores  D,  the  polarization  process  is  shut  off, 
limiting  the  total  nuclear  polarization  that  can  build  up.  While  not  shown  in  Fig.  2.5e, 
we  observe  a  dramatic  reduction  of  the  total  \D\,  not  just  Dz,  consistent  with  this 
qualitative  observation.  However,  because  we  have  not  observed  this  phenomenon  in 
any  physical  parameter  regimes  we  shall  not  study  it  further. 

Crossover 

For  many  choices  of  parameters,  we  find  both  trajectories  in  which  Dz  — »  0  and 
\DZ\  remains  large,  depending  on  initial  conditions,  as  shown  in  Fig.  2.6a.  Note  that 
when  we  add  a  small  amount  of  transverse  dephasing  to  these  trajectories,  as  shown 
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Figure  2.6:  a)  1000  trajectories  were  run  with  initial  conditions  chosen  from  the 
thermal  distribution  with  no  noise.  The  mean  value  of  \DZ\  is  shown  in  black,  and 
the  the  gray  region  enclosing  67%  of  the  trajectories.  A  single  trajectory  is  shown  in 
the  thin  red  line.  For  parameters,  see  Table  1.  These  parameters  are  not  represented 
in  the  phase  diagram  since  they  have  very  large  A+.  For  these  parameters,  many 
trajectories  are  attracted  near  D  =  0,  as  in  the  single  trajectory  shown,  for  extended 
periods  of  time,  b)  Trajectories  were  begun  from  identical  configurations  as  in  a,  this 
time  with  noise  added.  With  noise  included,  the  metastability  of  the  zero  state  is 
removed,  and  the  gray  region  is  now  bounded  away  from  zero. 
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in  Fig.  2.6b,  the  median  value  of  \DZ\  does  not  markedly  change,  but  there  are  no 
longer  trajectories  with  Dz  — >  0;  the  noise  apparently  disrupts  the  fragile  attraction 
toward  \DZ\  — »  0.  Simulations  performed  with  parameters  intended  to  approximate 
experiments  (Reilly  et  ah,  2008b;  Foletti  et  ah,  2009)  are  in  this  crossover  regime. 

Stability  of  Zero  States 

We  now  investigate  more  carefully  the  stability  of  the  zero  states.  Near  the  zero 
state  the  EOM  are  greatly  simplified  because  many  of  the  terms  in  arise  from 
perturbative  processes  involving  multiple  applications  of  D.  Keeping  only  the  terms 
linear  in  D  and  working  to  first  order  in  mo  we  can  write 

D+  =  (r0  +  -)S*ZD+  +  ( TQm0S;S+  -  iA0S*+)Dz  (2.9) 

Dz  =  — Re[(r0  +  iA_)D+S*_]  -  T0m0S±  ■  S*±DZ,  (2.10) 

where  we  have  introduced  the  variable  S*  =  Ylkddld^kd/^-  Because  dS/dt,dS* /dt  ~ 
O(D),  we  can  neglect  the  time  dependence  of  S  and  S*  in  the  EOM  for  D  near  the 
zero  state.  After  a  long  time  the  system  becomes  polarized  so  that  S*  <C  0,  this 
allows  us  to  adiabatically  eliminate  D+  to  obtain 


(Fo  +  AW)  \S*\ 

(2.11) 

4  =  0  +  0(D2) 

(2.12) 

This  linear  stability  analysis  gives  no  conclusion  about  the  stability  of  the  zeros 
states.  This  result  implies  that  within  this  model  the  stability  of  the  zero  state  is 
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only  determined  at  higher  order.  This  is  a  little  surprising  because  at  first  glance 
Eq.  2.10  appears  to  have  an  attractive  force  towards  Dz  =  0.  This  arises  from  the 
same  mechanism  described  by  Stopa  et  al.  (2010);  however,  a  more  careful  treatment 
reveals  that  this  effect  actually  cancels.  Our  simulations  indicate  that  the  nonlinear 
corrections  make  the  zero  state  repulsive  in  the  experimentally  relevant  parameter 
regimes.  When  we  include  the  nuclear  spin  noise  we  shall  show  analytically  that  the 
system  is  repelled  from  the  zero  states. 

2.4.2  Effect  of  Nuclear  Spin  Noise 

Unequal  Dots 

Our  results  that  zero  states  are  unstable  to  the  growth  of  large  difference  fields 
in  the  presence  of  asymmetry  in  the  size  of  the  dots  and  nuclear  noise  can  be  be 
understood  in  the  following  heuristic  picture  first  given  by  Gullans  et  al.  (2010).  We 
assume  the  nuclear  spins  have  equal  spin  flip  rates  on  the  two  dots,  which  is  borne 
out  by  the  analytical  and  numerical  calculations  presented  below.  Then  the  build¬ 
up  of  the  total  Overhauser  held  Sz  is  proportional  to  —  (gg  +  gr),  where  g^r)  are 
the  effective  hyperhne  interactions  on  the  left  (right)  dot  and  the  negative  sign  arise 
because  nuclear  spins  are  hipped  down  in  the  experimental  cycles.  Similarly  Dz  grows 
as  —  (ge  —  gr)  so  that  the  ratio 

Dz/Sz  >  {gg  —  gr)/{gg  +  gr)-  (2-13) 

In  this  section  we  demonstrate  a  similar  result  within  our  full  model.  We  assume 
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Figure  2.7:  a)  Asymptotic  value  of  \DZ/SZ\  as  a  function  of  dot  asymmetry  with 
parameters  chosen  as  in  the  location  marked  with  an  x  in  Fig.  2.4,  strongly  in  the 
instability  regime.  The  horizontal  access  corresponds  to  the  left  dot  decreasing  in  size 
from  right  to  left,  which,  by  our  simple  argument,  should  result  in  a  positive  ratio 
of  Dz/Sz.  Trajectories  which  show  the  opposite  sign  indicate  a  competition  with  the 
coherent  instability  mechansim.  For  each  value  of  dot  asymmetry  R,  we  initiahzed 
fifty  runs  in  a  single  initial  spin  configuration  chosen  from  the  thermal  distribution 
(with  Dz  =  —0.72  and  Sz  =  —1.57).  We  plot  the  asymptotic  value  of  Dz/Sz.  The 
runs  that  ended  with  Dz/Sz  greater  (less)  than  0  shown  are  shown  as  red  (blue) 
points.  The  circles  (crosses)  indicate  the  mean  value  of  the  red  (blue)  points,  with 
error  bars  showing  the  standard  deviations.  The  solid  and  dashed  lines  are  given 
by  Eq.  2.24  and  Eq.  2.13,  respectively,  b)  As  in  (a),  with  parameters  chosen  in  the 
location  marked  with  an  o  in  Fig.  2.4,  strongly  in  the  saturation  regime.  Here  the 
sign  of  the  ratio  Dz/Sz  follows  what  is  expected  from  the  natural  asymmetry. 
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homogeneous  coupling  and  work  in  the  high  field,  large  J,  limit  where  we  can  set 
Ao  =  A_  =  0  in  Prf.  The  local  noise  processes  included  in  Eq.  2.7  give  rise  to  a 
mean  decay  of  the  collective  nuclear  spin  variables  and  associated  fluctuations  Tt(r)  , 
for  L+(R+),  defined  by  (J7z(^)  J7d'(t>))n  =  2 D,^Sdd'^(t  —  t').  The  semiclassical  EOM 
for  the  nuclear  spins  reduce  to 

L+  —  ge r0  LZ(L+  —  R+)/ 2  —  //  L+  +  \/2r/  J- e,  (2-14) 

L2  =  -|r0(Li-R±-L±),  (2.15) 

and  similarly  for  R,  where  rj  is  defined  in  Eq.  2.8.  From  Eq.  2.14,  we  see  that  if  we 
start  in  a  zero  state,  T,i  will  produce  a  fluctuation  in  D±,  and  the  contribution  to  Lz 
of  the  form  —  ggT0L2±  results,  in  the  long  time  limit,  in  Lz  <C  —  1  and  similarly  for 
Rz.  Thus,  \LZ/LZ\  1  and  we  can  treat  Lz,  Rz  as  static  to  find  ( L2±)n ,  (R±)n  and 
(Lj_  •  Rj_)n,  which  allow  us  to  find  the  slow  evolution  of  Lz,  Rz. 

In  particular,  assuming  Lz,  Rz  are  constant  we  can  write  the  closed  set  of  equations 
for  L+  and  R+ 

1  L+ 

y  R+ 


\ 

/ 


2 


/ 


gt.Lz  —geLz 


\ 


(  r  \ 


-v 


y  grRz  grRz  J 

\! 2ij 


\R+  / 


<  lC 


1  T  t 

•r  f 
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Introducing  the  variables 


we  find 


(2.17) 


(2.18) 

(2.19) 


here  75  =  —  T0(g£Lz+grRz)/2  >  0.  We  can  use  this  solution  to  calculate  ( L2±)n  ,  {R2±)n, 
and  ( Lj_  ■  R±)n ■  For  example  to  lowest  order  in  1  /LZ)  1/RZ 


PC  = 


X 


*v/9 
(i  +  p )2 

9e  +  9rP2 
2  r] 


+  (9e  +  9r)p2  2p(ge  -  grp ) 
27s  7s 


(2.20) 


where  we  have  defined  p  =  giLz/grRz,  g  =  (ge+gr)/ 2  and  used  the  fact  that  D2  =  g^/g 
in  our  units. 

Inserting  this  solution  into  the  EOM  for  Dz,  Sz  gives  reduced  EOM  for  the  slow, 
noise- averaged  evolution  of  Dz  and  Sz.  After  some  straightforward  manipulations  we 
arrive  at 
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where  Sf  =  (gtLx  +  grRz )/2  =  -7s/r0  and 

(1  -  /7)(1  -  R2)  (1  -R)3  N 

(1  -  i?)(l  +  i?)2  -(1  +  i?)(l  -  i?2)  } 


(2.22) 


and  i?  =  (jr! 9 (■  After  rescaling  time  to 


T  = 


9i  V  gt  Qr 
9  \Sl(t')\2 


(2.23) 


this  becomes  a  purely  linear  system  characterized  by  the  matrix  E.  For  all  R  >  0,  this 
matrix  has  one  positive  and  one  negative  eigenvalue;  thus,  it  has  one  growing  mode 
and  one  decaying  mode.  In  the  long  time  limit,  both  Sz  and  Dz  will  be  proportional 
to  their  overlap  with  the  growing  mode.  Thus  Dz/Sz  approaches  a  constant,  which 
is  easily  found  from  E  as 


Dz  1  -  R2 

Sz  2 R  +  \/4:R2  +  (1  —  R)A ' 


(2.24) 


In  Fig.  2.7  we  compare  this  result  and  Eq.  2.13  to  the  full  numerics  including 
all  the  parameters.  The  horizontal  access  corresponds  to  the  left  dot  decreasing  in 
size  from  right  to  left,  since  Dz/Sz  ~  (ge  —  gr)/(ge  +  gr)  according  to  our  simple 
argument  we  expect  this  to  result  in  a  positive  ratio  of  Dz/ Sz.  In  Fig.  2.7a,  however, 
we  see  that  for  small  asymmetry  gr/ge  >  0.5,  many  trajectories  have  the  opposite 
sign  indicates  that  in  this  regime  the  coherent  instability  mechanism  (which  does  not 
prefer  either  sign)  competes  with  the  natural  asymmetry.  For  larger  asymmetries 
9r/ 9t  <  0.5  all  trajectories  are  seen  to  follow  the  direction  of  the  natural  asymmetry. 
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Fig.  2.7b  shows  the  same  simulations  performed  in  the  saturation  regime.  As  there 
is  no  coherent  instability  mechanism  competing  with  the  dot  asymmetry,  the  sign  of 
Dz  is  determined  by  the  asymmetry  in  all  but  the  most  symmetric  dots.  Dz/Sz  is  in 
good  agreement  with  the  simple  prediction  given  by  Eq.  2.13  and  Eq.  2.24. 

Identical  Dots 

For  identical  dots  the  arguments  given  in  the  previous  subsection  break  down; 
however,  we  shall  now  show  that  for  certain  parameters  there  still  exists  a  mechanism 
for  self-consistent  growth  of  \DZ\.  Growth  of  \DZ\  requires  nonzero  D j_.  For  inter¬ 
mediate  held  and  exchange,  the  Ao,_  contributions  to  Pd  become  comparable  to  the 
r0  term.  In  particular,  the  A0Dzz  term  acts  as  a  source  term  for  D±  (see  Eq.  2.25). 
Consequently,  for  weak  enough  noise  D±  will  only  be  appreciable  when  |Ao-D2/r0S'2| 
is  appreciable,  which  provides  a  self-consistency  condition  for  the  continued  growth 
of  Dz. 

These  properties  of  identical  dots  can  be  seen  analytically  in  the  following  limiting 
case:  we  assume  a  wave  function  where  the  coupling  takes  two  values,  g\  #2,  V  and 
that  initially  —  g2Sz  g{  \DZ\  g1,  S±  ~  1  and  D±  ~  Dz/Sz  <C  1.  We  denote  the 
total  angular  momentum  of  nuclear  spins  in  dot  d  with  coupling  constant  gk  by  Jkd 
and  assume  J^d  ~  ~  J|d  <C  J{d  so  that  the  majority  of  the  polarization  resides  in 

the  strongly  coupled  spins.  We  can  write  a  closed  set  of  equations  for  the  evolution 
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of  D  and  S 


D+  =  gliA_SzD+  -  g1iA0DzS+ 

+92  8  i  AqDz(J^  +  J2r)/ 2  —  g^  8iA_  D+(J 2t  +  J|r)/2, 

S+  =  -gii(A0  -  A -)DZD+  +  g28i  A0Dz(j£e  -  Jfr)/ 2 
—  g2S  i  A-D+(  J2e  —  J2r)/2, 

Dz  =  ,9l  Im  [A_D+S_]  -  g2  5  Im  [A_D+(  +  J2"  )/2] , 

^  =  -giY,Dl  -  g2Slm[A_D+(J-f  -  J~)/ 2], 
j^d.  =  ±92'IAqDz  J^d  =p  g2iA_D+J2d  —  gJ2d  +  fd , 
j“2d  =  ±92  Im[A -D+J~d\ , 

where  the  top  sign  is  for  d  —  £,  A_  =  A_  —  ir0,  8  =  g±  —  g2,  fd  is  a  gaussian, 
white  noise  process  derived  analogously  to  Td  such  that  ( fdfd)n  =  2 //<x2,  and  we 
have  neglected  to  write  the  noise  terms  in  the  EOM  for  D+  and  S+  because  we  have 
assumed  they  are  higher  order.  Furthermore,  we  can  neglect  all  terms  proportional 
to  g2D+Jfd  because  these  are  second  order.  This  leads  to  the  somewhat  simpler  set 
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of  equations 


D+  =  g1iA_SzD+  -  g1iA0DzS+ 

(2.25) 

+  g2  8  i  A0DZ  (Jj  +  J2r)/ 2, 

S+  =  -gii( A0  -  A -)DZD+ 

(2.26) 

+  g28iA0Dz(J+-J+)/2r 

j2d  =  T  g  2'i  A  ()DZ,  —  gJ2d  +  fd, 

(2.27) 

Dz  =  gi  Im[A_JD+S'_] 

(2.28) 

Sz  =  -giT0Dl, 

(2.29) 

These  equations  can  be  solved  perturbatively  in  l/Sz,l/Dz  by  the  same  method  as 
in  the  previous  section.  The  only  difference  in  the  structure  of  the  two  problems  is 
that  in  this  case  the  source  terms  for  D+  and  S+  are  proportional  to  J+d  instead  of 
white  noise;  as  a  result  we  have  to  take  into  account  the  coherent  evolution  of  the 
source  term.  We  can  expand  the  resulting  EOM  for  Dz  in  g\Dzjg2Sz  to  find  the 
noise-averaged  equation 

D«  =  -9lr02iV2(  filv) 

W  +  (2'30) 
q  +  A2_  A I  Al/ 

from  which  we  see  that  the  sign  of  Tg  +  A2_  —  AoA_  determines  whether  or  not  there 
is  continued  growth  of  Dz.  Note  that  the  perturbation  theory  breaks  down  as  g2  — >  0. 
This  reflects  the  importance  of  including  the  coherent  evolution  of  J+d  in  solving  for 
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Figure  2.8:  Phase  diagram  as  in  Fig.  2.4  except  with  noise  added.  The  phase  diagram 
is  nearly  identical.  See  Table  I  for  parameters. 


the  dynamics.  Without  5%  we  would  have  found  Dz  =  0.  This  phase  boundary  is 


shown  as  the  dashed  line  in  Fig.  2.3.  In  Fig.  2.8  we  show  the  phase  diagram  as  a 
function  of  cycle  frequency  and  inverse  magnetic  field,  where  we  see  qualitatively  the 
same  behavior  as  Fig.  2.4. 

2.5  Relevance  to  Other  Central  Spin  Systems 

Although  this  work  has  focused  on  lateral  double  quantum  dots  in  GaAs,  the 
methods,  and  some  of  the  results,  can  be  applied  to  vertical  double  dots  (Takahashi 
et  ah,  2011),  InAs  quantum  dots  (Sun  et  al.,  2012;  Hogele  et  ah,  2012),  silicon  based 
quantum  dots  (Maune  et  ah,  2012),  and  NV-centers  in  diamond  (Childress  et  ah, 
2006).  A  few  important  differences  for  these  other  central  spin  systems  are  that  the 
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sign  of  the  electron  ^-factor  may  be  positive  (compared  to  its  negative  sign  in  GaAs) 
and  the  spin-orbit  coupling  can  be  much  larger  in  other  systems  than  it  is  in  GaAs 
(Stepanenko  et  ah,  2012).  The  results  presented  in  the  paper  are  not  dependent  on 
the  sign  of  the  g-faetor.  Changing  the  sign  would  reverse  the  direction  of  the  nuclear 
polarization  from  negative  to  positive,  but  all  of  our  analysis  would  carry  through 
essentially  unchanged.  The  competition  between  spin-orbit  coupling  and  DNP  is 
more  dramatic  and  can  have  a  qualitative  effect  on  the  polarization  dynamics  for 
large  spin-orbit  coupling  (Rudner  and  Levitov,  2010). 

2.6  Conclusions 

We  have  shown  that  dynamic  nuclear  polarization  experiments  in  double  quantum 
dots  give  rise  to  a  rich  set  of  phenomena.  We  find  that  after  many  thousands  of  nuclear 
spin  pumping  cycles,  corresponding  to  experimental  timescales  of  several  hundred 
microseconds,  the  total  nuclear  polarization  is  driven  to  10  —  30%  of  full  polarization. 
The  polarization  is  aligned  opposite  the  magnetic  field  as  opposed  to  the  thermal 
polarization.  In  addition  to  this  large  polarization,  we  find  the  competition  between 
polarization,  noise  processes  and  coherent  evolution  mediated  by  the  electrons  allows 
one  to  carefully  control  the  final  nuclear  spin  state  in  the  two  dots.  We  have  developed 
detailed  numerical  and  analytical  methods  to  theoretically  describe  such  dynamics; 
however,  our  analysis  is  semiclassical  and  leaves  out  effects  such  as  spin-orbit  coupling 
and  a  full  description  of  the  nuclear  dipole-dipole  interactions  (which  we  approximate 
as  nuclear  spin  noise),  both  of  which  may  be  important  for  a  complete  understanding 
of  the  experiments. 
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The  main  implication  of  the  paper  for  DNP  experiments  in  double  dots  is  that  the 
nuclear  spin  dynamics  are  dominated  by  either  rapid  saturation  of  polarization  or  an 
instability  to  the  growth  of  large  difference  fields.  These  results  are  consistent  with 
the  experimental  observations  reported  in  Petta  et  al.  (2008),  Foletti  et  al.  (2009)  and 
Barthcl  et  al.  (2012);  however,  we  see  evidence  that  the  dynamics  are  much  richer  as 
the  experiments  have  not  resolved  whether  or  not  the  instability  to  large  difference 
fields  results  from  dot  asymmetry  or  coherent  electron-nuclear  interactions.  These  two 
cases  could  be  experimentally  distinguished  by  measuring  the  sign  of  Dz  in  a  given 
double  dot.  Furthermore,  we  showed  that  the  zero  states  may  be  experimentally 
observable  as  metastable  states  in  certain  parameter  regimes,  indicating  that  there  is 
still  much  to  explore  in  the  polarization  dynamics  of  double  quantum  dots. 
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Chapter  3 

Nanoplasmonic  Lattices  for 
Ultracold  Atoms 

3.1  Introduction 

Coherent  optical  fields  provide  a  powerful  tool  for  manipulating  ultracold  atoms 
(Bloch  et  ah,  2012;  Grimm  et  ah,  2000).  However,  diffraction  sets  a  fundamental  limit 
for  the  length-scale  of  such  manipulations,  given  by  the  wavelength  of  light  (Hecht, 
1998).  In  particular,  the  large  period  of  optical  lattices  determines  the  energy  scale  of 
the  associated  many-body  atomic  states  (Buluta  and  Nori,  2009;  Yi  et  ah,  2008;  Leung 
et  ah,  2012;  Lewenstein  et  ah,  2012).  The  resulting  scaling  can  be  best  understood 
by  noting  that  in  the  first  Bloch  band  the  maximum  atomic  momentum  ~  l/£,  where 
t  is  the  lattice  spacing.  This  sets  the  maximum  kinetic  energy  to  h2 /ml2  (Jaksch 
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et  al.,  1998).  For  conventional  optical  lattices  the  lattice  spacing  is  set  by  half  the 
wavelength  of  the  trapping  light  ~  500  nm;  this  yields  corresponding  tunneling  rates 
of  up  to  a  few  tens  of  kHz.  Additionally,  for  atoms  in  their  electronic  ground  states 
interactions  are  restricted  to  short  range. 

Recent  experimental  (Stehle  et  ah,  2011)  and  theoretical  (Murphy  and  Hau,  2009; 
Chang  et  ah,  2009)  work  has  demonstrated  that  integrating  plasmonic  systems  with 
cold  atoms  represents  a  promising  approach  to  achieving  subwavelength  control  of 
atoms.  In  particular,  the  experiments  of  Ref.  (Stehle  et  ah,  2011)  showed  that  ul¬ 
tracold  atoms  can  be  used  to  probe  the  near  fields  of  plasmonic  structures,  paving 
the  way  to  eventually  trap  atoms  above  such  structures.  In  this  chapter  we  propose 
and  analyze  a  novel  approach  to  the  realization  of  high-density  optical  lattices  us¬ 
ing  the  optical  potential  formed  from  the  near  held  scattering  of  light  by  an  array 
of  plasmonic  nanoparticles.  By  bringing  atom  trapping  into  the  subwavelength  and 
nanoscale  regime  we  show  that  the  intrinsic  scales  of  tunneling  and  onsite  interaction 
for  the  Hubbard  model  can  be  increased  by  several  orders  of  magnitude  compared  to 
conventional  optical  lattices.  In  addition,  subwavelength  confinement  of  the  atoms 
results  in  strong  radiative  interactions  with  the  plasmonic  modes  of  the  nanoparticles 
(de  Leon  et  ah,  2012).  The  coupled  atom-plasmon  system  can  be  considered  as  a 
scalable  cavity  array  that  results  in  strong,  long  range  spin-spin  interactions  between 
the  atoms  with  both  dissipative  and  coherent  contributions  (Cirac  et  ah,  1997;  Kim¬ 
ble,  2008).  Such  a  system  can  be  used  for  entanglement  of  remote  atoms  as  well  as 
for  novel  realizations  of  coherent  and  dissipative  many-body  systems. 
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Figure  3.1:  a)  Illustration  of  the  relevant  physics  in  the  plasmonic  lattice,  b)  Il¬ 
lustration  of  how  to  engineer  a  blue-detuned  optical  dipole  trap  by  driving  on  the 
blue  side  of  the  plasmon  resonance,  c)  Atomic  potential  for  Rb  including  van  der 
Waals  (vdw)  for  trapping  above  a  single  silver  nanoshell.  Dotted  line  shows  how  to 
weaken  the  trap  by  applying  circularly  polarized  light  perpendicular  to  the  trapping 
light.  (Inset)  Real  (dashed)  and  imaginary  (solid)  part  of  the  dipole  polarizability 
for  a  sphere  and  the  nanoshell  with  a  15  run  radius  and  13.85  nm  Si02  core,  d)  y-z 
contours  of  atomic  potential  in  MHz  for  a  line  of  nine  spheres  in  the  center  of  a  45x45 
square  lattice  with  a  60  nm  lattice  spacing,  black  regions  are  where  the  potential  is 
negative  due  to  vdw,  spheres  are  shown  in  white.  The  nanoshells  are  silver  with  a  15 
nm  radius  and  13.65  nm  Si02  core,  the  trapping  light  is  red  detuned  (704  nm)  wrt  to 
the  plasmon  resonance  (682  nm)  and  applied  from  above  with  rotating  x-y  polarized 
light. 
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3.2  Atom  Trapping  Above  a  Single  Metallic  Nanopar¬ 
ticle 

To  illustrate  our  approach  we  first  consider  a  single  metallic  nanosphere  in  vac¬ 
uum  illuminated  by  a  plane  wave.  For  spheres  small  compared  to  a  wavelength  the 
dominant  contribution  to  the  scattered  field  is  the  dipole  term,  where  the  induced 
dipole  moment  is  given  by  p  =  a(uj)E0  with 

a(u)  =  4vre0a3^|— ]-  (3.1) 

£{c a)  +  2 


where  a  is  the  radius  of  the  sphere  and  £  is  the  permitivity  (Jackson,  1999).  The  total 


electric  field  is 


E  —  E{\ 


a(u)  3 (r  ■  E0)r  —  E0 


Near  £(ujsp)  =  —2  there  is  a  plasmon  resonance  and  the  scattered  field  can  be 
engineered  to  create  an  optical  dipole  trap  as  depicted  in  Fig.  3.1b.  Specifically, 
when  the  applied  field  is  linearly  polarized  on  the  blue  side  of  the  plasmon  reso¬ 
nance  then  the  induced  dipole  will  be  ~  n  out  phase  with  the  incident  field,  lead¬ 
ing  to  two  intensity  minima  along  the  polarization  direction  at  the  positions  zf  = 
±2a3u2p/ (u2  —  (j o2p),  where  we  took  a  Lorentzian  polarizability  near  the  resonance 
a(uj)  =  4:7Te0a3uj2p/ (lo2p  —  u>2  —  icon),  with  n  the  linewidth.  For  red  detuned,  circu¬ 
larly  polarized  light,  there  will  be  two  minima  along  the  propagation  axis.  An  atom 
can  be  trapped  in  these  intensity  minima  via  optical  dipole  forces  (Grimm  et  ah, 
2000).  The  trapping  potential  is  given  by  htt2/5,  where  12  =  pi0  ■  E/h  is  the  Rabi 
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frequency,  /i0  is  the  atomic  dipole  moment,  and  6  =  ua  —  u>  is  the  detuning  between 
the  atom  and  laser.  Expanding  near  the  trap  minima  gives  the  trapping  frequency 
uif  =  9  j^2  Re  (a)2/  |or|2  ~  HQq/5  ma2. 

The  trap  depth  can  be  controlled  by  applying  a  second  held  with  the  opposite 
polarization,  as  illustrated  in  Fig.  3.1c.  Using  this  method,  the  atoms  can  be  loaded 
into  the  near  held  traps  by  starting  with  a  cold,  dense  gas  of  atoms  in  a  large  trap 
and  then  adiabatically  turning  on  the  near  held  traps. 

We  now  address  several  practical  considerations.  First,  for  alkali  atoms  there  is 
a  large  disparity  between  the  natural  plasmon  resonance  and  the  atomic  trapping 
transitions.  For  a  solid  silver  sphere  the  plasmon  resonance  occurs  near  350  nm 
(Johnson  and  Christy,  1972),  compared  to  780  nm  for  the  D2  line  in  Rb.  However, 
the  plasmon  resonance  is  easily  tuned  by  changing  the  geometry.  Adding  an  inert 
core,  such  as  SiCU,  will  shift  the  plasmon  resonance  into  the  red  (Bohren  and  Huffman, 
1983),  as  illustrated  in  the  inset  to  Fig.  3.1c. 

There  will  also  be  significant  surface  interactions.  In  Appendix  B.l  we  calculate 
the  Van  der  Waals  (vdw)  interaction  between  a  single  metal  sphere  and  an  atom. 
These  vdw  forces  can  be  overcome  with  relatively  modest  laser  power  because  of 
the  sphere’s  plasmonic  enhancement  (Murphy  and  Hau,  2009;  Chang  et  ah,  2009). 
There  are  two  dominant  sources  of  heating  and  decoherence  arising  from  incoher¬ 
ent  transitions  induced  by  the  trapping  laser  or  thermal  magnetic  held  noise  in  the 
nanoparticle.  The  first  effect  scales  as  yf 12/J2,  where  7  is  the  atomic  linewidth, 
and  is  suppressed  at  large  detuning.  To  estimate  the  effect  of  magnetic  held  noise 
we  approximate  the  nanoshell  as  a  current  loop  of  radius  and  height  a,  thickness 
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t,  and  resistivity  p.  Then  the  incoherent  transition  rate  between  hyperhne  states  is 
~  (gFPoUB)2kBT(a4t/r5)/h2pr,  where  r  is  the  distance  of  the  atom  to  the  sphere  cen¬ 
ter,  gp  is  the  hyperhne  g-factor,  ps  is  the  Bohr  magneton,  and  T  is  the  temperature 
(Henkel  et  ah,  1999). 

3.3  Atom  Trapping  Above  a  Lattice  of  Nanoparti¬ 
cles 

Figures  3. led  show  the  atomic  trapping  potential  for  a  single  sphere  and  an  array, 
respectively.  We  numerically  obtained  the  trapping  potential  in  Fig.  3.1c  using  Mie 
theory  and  the  vdw  potential  was  obtained  using  the  methods  in  Ref.  (Reid  et  ah, 
2009).  To  solve  for  the  trapping  potential  in  the  array  in  Fig.  3.  Id  we  approximated 
the  scattered  held  from  each  nanoshell  by  a  dipole  and  solved  for  the  total  held 
sclf-consistently.  Using  the  parameters  in  Fig.  3.1c  for  trapping  8'Rb  above  a  silver 
nanoshell  at  room  temperature  with  Ho  =  25  GHz  (corresponding  to  ~  108/sat,  where 
Gat  ~  1-7  mW/cm2)  and  5  =  25  THz,  we  estimate  a  trap  depth  of  ~  25  MHz  and 
a  trapping  frequency  of  ~  5  MHz.  Both  the  magnetic  held  noise  and  laser  detuning 
limit  the  decoherence  rate  to  ~  10  Hz  and  the  heating  rate  to  ~  1  Hz,  meaning  that 
the  atom  can  be  trapped  for  ~  1  second. 

The  controlled  patterning  of  arrays  of  metallic  nanoparticles  can  be  done  litho¬ 
graphically  in  a  top-down  approach  or  through  the  controlled  self-assembly  of  metallic 
nanoparticles  in  a  bottom-up  approach  (Nagpal  et  al.,  2009;  Lindquist  et  ah,  2012; 
Fan  et  ah,  2010;  Grzelczak  et  ah,  2010).  In  any  nanofabricated  system  one  must 
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contend  with  disorder;  the  relevant  disorder  in  this  system  occurs  in  the  particle 
positioning  and  particle  formation.  In  lithographic  approaches  one  can  control  the 
particle  formation  at  the  level  of  1-2  nm  (Lindquist  et  ah,  2012).  In  bottom- up,  self- 
assembly  approaches  it  is  possible  to  create  large  regions  of  well  ordered  crystal  with 
a  finite  density  of  point  and  line  defects,  much  like  a  conventional  solid  (Grzelczak 
et  ah,  2010).  Due  to  the  local  nature  of  the  traps  the  disorder  in  the  particle  posi¬ 
tioning  will  not  affect  the  trapping.  Errors  in  the  particle  formation  can  influence  the 
trap  by  shifting  the  plasmon  resonance  and  the  held  enhancement  of  each  particle. 
To  achieve  consistent  traps  the  fractional  error  in  the  plasmon  resonance  should  be 
smaller  than  its  inverse  quality  factor  Q  =  cosp/k,  which  for  silver(gold)  nanospheres 
goes  up  to  80(20)  (Johnson  and  Christy,  1972;  Hartland,  2011).  Currently,  metallic 
nanoshells  can  be  made  with  a  fractional  error  in  the  radius  of  less  than  5%,  which 
is  comparable  to  the  inverse  of  Q  (Rycenga  et  al.,  2011). 

3.4  Hubbard  Models  in  Nanoscale  Lattices 

As  a  first  example  application  of  this  system  we  consider  a  realization  of  the 
single-band  Hubbard  model  in  the  novel  regime  of  large  atomic  density  (Bloch  et  ah, 
2012).  As  an  example,  Fig.  3.  Id  shows  that  a  well  defined  lattice  potential  can  be 
achieved  with  a  period  of  60  nm,  which  is  within  current  fabrication  limits.  Figure 
3.2a  illustrates  the  scaling  for  the  maximum  tunneling  in  the  lowest  band  and  the 
corresponding  on-site  interaction  Uq.  (Jaksch  et  al.,  1998).  In  Appendix  B.3  we  show 
that  the  tunneling  rate  can  also  be  tuned  through  appropriate  polarization  control. 

These  nanoscale  traps  reach  a  regime  of  atomic  confinement  where  the  ground 
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state  uncertainty  becomes  comparable  to  the  free  space  scattering  length.  For  two 
atoms  in  a  3D  isotropic  trap  the  two-body  scattering  problem  can  be  solved  exactly, 
leading  to  an  effective  scattering  length  aes(coT )  which  depends  on  the  confinement 
energy  (Busch  et  al.,  1998;  Bolda  et  ah,  2002).  The  inset  of  Figure  3.2  shows  that  a 
resonance  emerges  in  the  effective  scattering  length  as  a  function  of  trap  frequency. 
We  show  how  this  is  calculated  in  Appendix  B.4. 

Disorder  in  the  lattice  will  also  effect  the  Hubbard  model.  The  dominant  effect 
arises  from  shifts  in  the  local  atomic  potential  at  each  sphere  as  the  plasmonic  en¬ 
hancement  factor  changes  from  site  to  site.  From  Eq.  2  one  can  show  that  the  rms  of 
the  disorder  potential  is  given  by  U^is  ~  zf/a9Q2)r]/ujSp ,  where  rj  is  the  rms  error 
in  the  plasmon  resonance.  If  we  take  rj /ojsp  ~  5%,  then  for  a  wide  range  of  param¬ 
eters,  including  those  in  Fig.  3.  Id,  we  find  that  UdiS  can  be  made  smaller  than,  or 
comparable  to,  the  maximum  tunneling.  In  addition,  since  the  disorder  is  static  one 
can  reduce  it  using  the  techniques  described  in  Ref.  (Pichler  et  ah,  2012).  The  effect 
of  disorder  on  the  single-particle  physics  is  well  understood  (Lagendijk  et  ah,  2009); 
moreover,  the  interplay  between  interactions  and  disorder  in  the  Hubbard  model,  as 
studied  in  Ref.  (Bclitz  and  Kirkpatrick,  1994;  Basko  et  ah,  2006;  Byczuk  et  ah,  2005; 
Fallani  et  ah,  2007),  is  an  interesting  new  regime  which  can  be  explored  in  the  present 
system. 
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Figure  3.2:  Shows  the  scaling  of  the  maximum  tunneling  in  the  lowest  band,  and 
the  corresponding  on-site  interaction.  Calculated  using  the  Wannier  functions  for  a 
sinusoidal  potential.  (Inset)  Energy  dependent  scattering  length  for  two  87Rb  atoms 
on  a  single  site  as  a  function  of  the  trap  frequency. 

3.5  Plasmon  Mediated  Interactions  and  Entangle¬ 
ment  in  the  Nanolattice 

We  now  consider  long  range  interactions  within  the  plasmonic  lattice,  associated 
with  the  strong  radiative  coupling  between  the  atoms  and  spheres  (Genov  et  ah, 
2011).  This  can  be  viewed  as  a  strongly  coupled  cavity  QED  system.  The  coupling 
between  the  atoms  and  the  near  field  of  the  sphere  is  given  by  g  ~  /i0d0/e0r3  where 
do  =  \Jhu)spa( 0)/2  is  the  quantized  dipole  moment  of  the  sphere  (de  Vries  et  ah, 
1998).  Since  the  plasmons  are  overdamped  the  relevant  coupling  is  given  by  the 
Purcell  factor  P  =  g2 /k 7.  The  plasmon  linewidth  k  has  contributions  from  radiative 
and  ohmic  losses.  The  radiative  damping  rate  is  k3do/3ne0h  ~  k3a3cusp.  Large  spheres 
are  radiatively  broadened  and,  in  this  case,  P  ~  ( kr)~ 6,  while  for  small  spheres 
P  ~  Qa3/k3r6.  In  both  limits,  when  r  <C  A/27T  ~  100  nm  the  atoms  enter  the 
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Figure  3.3:  a)  Shows  the  cavity  QED  figure  of  merit  g2 /k 7  with  changing  system  size 
assuming  the  atom  is  trapped  at  a  distance  of  twice  the  sphere  radius.  We  show  the 
scaling  for  both  silver  and  gold  nanoshells  with  a  Q  of  80  and  20,  respectively.  (Inset) 
Single  atom  trapped  above  a  nanosphere  acts  as  cavity  QED  system  with  atomic 
and  cavity  losses  7  and  k,  respectively,  and  a  coherent  coupling  g.  b)  Fidelity  for 
generating  a  ground  state  singlet  state  between  two  atoms  on  the  lattice  with  their 
separation  after  optimization.  The  entanglement  is  generated  through  interaction 
with  the  collective  plasmon  modes,  where  we  took  the  metal  losses  of  bulk  silver. 
(Inset)  Scalable  cavity  QED  array  of  atoms  and  plasmons. 


strong  coupling  regime  P  1,  see  Figure  3.3a.  Note  that  there  are  also  multipolar 
corrections  to  the  Purcell  factor,  but  in  Appendix  B.l-2  we  show  these  scale  as  Im((e— 
l)/(e  +  l))a5/r5  ~  10-4  for  silver. 

For  a  lattice  of  nanospheres,  intersphere  coupling  is  also  present  and  leads  to 
delocalized  plasmon  modes  in  the  lattice  (Quinten  et  al.,  1998;  Krenn  et  al.,  1999). 
We  calculate  the  interaction  of  two  atoms  through  these  modes  in  a  ID  chain  of 
nanospheres.  For  each  sphere  in  the  chain  we  can  write  the  self-consistent  equation 
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for  their  dipole  moments  as  (Park  and  Stroud,  2004) 

Pn  4“  NnmPrn)  (3-3) 

where  pn  is  the  induced  dipole  moment  of  the  nth  nanopoarticle,  En  is  the  incident 
held,  and  Nnm  is  the  3x3  matrix  that  gives  the  dipole  held  at  site  n  due  to  the 
dipole  at  site  m.  In  ID  two  sets  of  transverse  modes  where  the  dipoles  are  oriented 
perpendicular  to  the  chain  and  one  set  of  longitudinal  modes  for  parallel  orienta¬ 
tion.  Defining  pq  to  be  the  qth.  eigenvector  of  Nnm  with  eigenvalue  Dq,  then  the 
effective  polarizability  of  the  qth  mode  is  af 1  =  or1  —  Dq ,  i.e.  pq  =  aqEq.  For  a 
Lorentzian  polarizability  the  real  part  of  Dq  gives  the  shift  in  the  resonance  frequency 
of  the  qth  mode  and  the  imaginary  part  gives  the  change  in  the  linewidth.  Nnm  is 
diagonalized  by  Fourier  transform  and  if  we  neglect  all  but  nearest  neighbor  terms 
Dq  =  2 jVgj  cos q  —  ik3/6ne0,  where  =  Re(-/V0i). 

Let  us  consider  atoms  trapped  above  the  ID  array  of  spheres.  The  plasmonic 
modes  can  be  adiabatically  eliminated  using  standard  methods  in  quantum  optics 
(Gross  and  Haroche,  1982).  For  two-level  atoms  polarized  parallel  to  the  ID  chain 
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the  atomic  density  matrix  evolution  is 


P  = 


lUat 
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(3.4) 

(3.5) 

(3.6) 


where  z  is  the  position  of  the  atoms  above  the  sphere  and  q*  =  q*  +  i  q*  is  the 
resonant  wavevector  such  that  a~}{ua)  =  0.  The  first  line  in  Eq.  3.4  describes  the 
coherent  evolution  and  the  second  line  describes  the  collective  dissipation.  Here  we 
have  neglected  the  contribution  to  the  interaction  from  free-space  radiative  modes. 

The  coherent  and  dissipative  contributions  to  Eq.  3.4  are  equally  strong  when  the 
atom  and  plasmon  are  near  resonant.  Working  far  off  resonance,  however,  results  in 
purely  coherent  dynamics,  which  can  be  used  to  implement  long-range  interacting 
spin  models  including  frustration  (Strack  and  Sachdev,  2011;  Gopalakrishnan  et  ah, 
2009;  Gardner  et  al.,  2010).  Alternatively,  the  collective  dissipative  dynamics  can  be 
used  to  prepare  correlated  atomic  states  (Verstraete  et  ah,  2009).  As  an  example,  we 
now  show  how  to  directly  prepare  a  ground  state  singlet  between  two  atoms  separated 
by  large  distances  on  the  lattice.  We  take  two  ground  states  | g)  and  |s)  and  an  excited 
state  |e)  which  is  coupled  to  | g)  via  an  external  field  and  only  decays  via  the  plasmons 
back  to  | g)  (see  inset  to  Fig.  3.3a).  An  external  microwave  field  mixes  the  two  ground 
states.  To  prepare  the  singlet  state  \S)  =  \gs)  —  \sg)  we  use  a  similar  approach  to 
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Ref.  (Kastoryano  et  ah,  2011)  whereby  the  singlet  state  is  engineered  to  be  the  steady 
state  of  a  driven,  dissipative  evolution.  We  take  a  separation  n  such  that  cos  q*n  =  1 
and 

P  =  -AovVWr  +  <kT]p  -  ^n{V[afe]  +  V[a92e])p  (3.7) 

where  V[c\p  =  l/2{cdc,  p}  —  cpC  and  8^n  =  700  —  7no  ~  7oo  ( i3/a3)n/Q  for  n  <C  Q. 
The  dynamics  can  be  mapped  to  a  cavity  QED  system  by  identifying  70 n  with  the 
collective  decay  g2 / k  and  8^n  with  the  free  space  decay  7.  The  two  excited  states  | eg) 
and  | ge)  split  into  a  superradiant  state  | eg)  +  | ge)  and  a  subradiant  state  | eg)  —  \ge) 
with  decay  rates  270n  +  8^n  and  8jn,  respectively. 

The  singlet  preparation  proceeds  as  follows.  First,  we  selectively  excite  the  sub¬ 
radiant  transition  \gg)  to  | ge)  —  \ eg)  by  driving  with  a  weak  external  laser  held 
Q  ~  700,  which  we  take  to  have  a  71  phase  difference  on  the  two  atoms.  Second, 

in  order  to  make  the  singlet  state  a  unique  steady  state,  we  apply  a  global  microwave 
held  to  mix  the  triplet  ground  states  without  affecting  the  singlet  state.  In  the  result¬ 
ing  dynamics,  the  pumping  rate  into  the  singlet  state  is  Ut2 /8^n)  while  the  pumping 
rate  back  into  the  triplets  is  fl2/7oo  (see  Appendix  B.5  for  more  details).  The  steady 
state  of  this  process  gives  the  singlet  state  with  hdclity  F  =  (/S'!  p  |  -S')  1  -  1/P' 

where  P'  =  7oo/^7n-  Fig.  3.3b  shows  the  hdelity  for  two  atoms  with  variable  separa¬ 
tion  obtained  from  numerical  simulation  of  Eq.  3.4. 

To  measure  the  correlations  in  this  system,  an  all  optical  approach  could  be  re¬ 
alized  by  making  the  nanoparticle  array  in  the  near  held  of  a  solid  immersion  lens 
(SIL),  which  enhances  the  resolution  beyond  the  diffraction  limit  by  a  factor  of  n, 
the  index  of  refraction  of  the  SIL  (Wu  et  ah,  1999).  Combining  a  SIL  with  e.g.  super 
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resolution  microscopy  techniques  would  allow  one  to  reach  the  requisite  resolution  of 
~50  nm  at  optical  wavelengths  (Huang  et  ah,  2009). 

3.6  Conclusions 

Our  analysis  shows  that  combining  cold  atom  techniques  with  nanoscale  plas- 
monics  reaches  new  regimes  in  controlling  both  the  collective  motion  of  atoms  and 
atom-photon  interactions.  Combining  excellent  quantum  control  of  isolated  atoms 
with  nanoscale  localization,  may  open  up  exciting  new  possibilities  for  quantum  con¬ 
trol  of  ultracold  atoms. 
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Chapter  4 

Single  Photon  Nonlinear  Optics 
with  Graphene  Plasmons 

4.1  Introduction 

Nonlinear  optical  processes  find  ubiquitous  use  in  modern  scientific  and  technolog¬ 
ical  applications,  facilitating  diverse  phenomena  like  optical  modulation  and  switch¬ 
ing,  spectroscopy,  and  frequency  conversion  (Boyd,  2003).  A  long-standing  goal  has 
been  to  realize  nonlinear  effects  at  progressively  lower  powers,  which  is  difficult  given 
the  small  nonlinear  coefficients  of  bulk  optical  materials.  The  ultimate  limit  is  that 
of  single-photon  nonlinear  optics,  where  individual  photons  strongly  interact  with 
each  other.  Realization  of  such  nonlinear  processes  would  not  only  facilitate  peak 
performance  of  classical  nonlinear  devices,  but  also  create  a  unique  resource  for  im¬ 
plementation  of  quantum  networks  (Kimble,  2008)  and  other  applications  that  rely 
on  the  generation  and  manipulation  of  non-classical  light. 
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One  approach  to  reach  the  quantum  regime  involves  coupling  light  to  individual 
quantum  emitters  (Duan  and  Monroe,  2008;  Kimble,  2008),  in  order  to  take  advantage 
of  their  intrinsically  nonlinear  electronic  spectrum.  While  a  number  of  remarkable 
phenomena  have  been  demonstrated  using  these  systems  (Haroche,  2013),  their  re¬ 
alization  remains  a  challenging  task.  Specifically,  in  contrast  to  conventional  bulk 
nonlinear  systems,  coherent  single  quantum  emitters  are  generally  unable  to  oper¬ 
ate  under  ambient  conditions,  suffer  from  relatively  slow  operating  speeds,  and  have 
limited  tunability  of  their  properties. 

Motivated  by  these  considerations,  recently  there  has  been  renewed  interest  in  bulk 
nonlinear  systems  that  can  reach  the  quantum  regime  (Matsuda  et  ah,  2009;  Mabuchi, 
2011;  Ferretti  and  Geraee,  2012).  In  particular,  recent  experiments  demonstrated  re¬ 
alization  of  a  quantum  nonlinear  medium,  featuring  single  photon  blockade  (Peyronel 
et  ah,  2012)  and  conditional  nonlinear  two-photon  phase  shifts  (Peyronel  et  ah,  2013), 
based  on  strongly  interacting  ultracold  atoms.  The  essence  of  this  approach  is  that 
the  probability  for  two  photons  to  interact  can  become  substantial  if  the  photons  are 
confined  to  a  sufficiently  small  mode  volume  of  the  nonlinear  medium  for  sufficiently 
long  times.  Motivated  by  these  recent  developments,  in  this  Letter  we  explore  the 
potential  for  using  nanoscale  surface  plasmon  excitations  in  graphene  (Mikhailov  and 
Ziegler,  2007;  Jablan  et  ah,  2009)  for  quantum  nonlinear  optics.  In  particular,  recent 
theoretical  (Mikhailov  and  Ziegler,  2007;  Jablan  et  ah,  2009;  Koppens  et  ah,  2011) 
and  experimental  (Fei  et  ah,  2012;  Chen  et  ah,  2012)  results  indicate  that  graphene 
plasmons  can  be  confined  to  volumes  millions  of  times  smaller  than  the  diffraction 
limit.  We  show  that  under  realistic  conditions,  this  field  confinement  enables  deter- 
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ministic  interaction  between  two  plasmons  (i.e.,  photons)  over  picosecond  time  scales, 
which  is  much  shorter  than  the  anticipated  plasmon  lifetime  (Principi  et  ah,  2013). 
We  show  how  one  can  take  advantage  of  this  interaction  to  realize  a  single  photon 
switch  and  produce  non-classical  light. 

4.2  Graphene  Plasmonics 

Graphene,  a  single  atomic  layer  of  carbon  atoms,  has  attracted  tremendous  in¬ 
terest  for  its  unique  electronic,  mechanical,  and  quantum  transport  properties  (Geim 
and  Novoselov,  2007;  Castro  Neto  et  al.,  2009).  More  recently,  its  optical  response 
has  also  been  explored.  For  example,  it  has  been  demonstrated  that  the  graphene 
band  structure  yields  a  constant  attenuation  rate  of  light  through  a  single  layer  of 
7tq;  ta  2.3%  when  the  graphene  is  in  its  intrinsic  (undoped)  state,  where  a  ~  1/137 
is  the  fine-structure  constant  (Nair  et  ah,  2008).  The  band  structure  also  produces 
remarkable  properties  for  guided  electromagnetic  surface  waves  in  the  form  of  surface 
plasmons  (SPs)  (Wunsch  et  ah,  2006),  as  we  now  describe. 

Through  electrostatic  gating,  it  is  possible  to  introduce  a  net  carrier  concentration, 
which  shifts  the  Fermi  energy  Hujf  away  from  the  Dirac  point  to  a  non-zero  value. 
The  in-plane  conductivity  of  graphene  is  well-approximated  by  the  expression  <j(uj)  a 
at  frequencies  below  twice  the  Fermi  frequency  u  <  2up  (Falkovsky,  2008), 
which  describes  a  Drude-like  response  of  electrons  within  a  single  band.  In  realistic 
systems  the  conductivity  will  also  have  a  small  term  7  describing  dissipation  due 
impurity  or  phonon-mediated  scattering.  There  are  two  limits  on  the  existence  of 
low-loss  SP  modes  in  graphene.  First,  at  frequencies  oj  >  2c up,  graphene  suffers  from 
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strong  inter-band  absorption  (e.g.,  giving  rise  to  its  attenuation  of  2.3%)  (Jablan 
et  al.,  2009;  Koppens  et  al.,  2011).  Second  for  frequencies  above  the  the  optical  phonon 
frequency  huop  ~  0.2  eV,  there  is  additional  loss  due  to  scattering  into  optical  phonons 
(Jablan  et  ah,  2009;  Yan  et  ah,  2013).  With  this  in  mind,  we  focus  on  the  regime 
where  the  frequencies  fall  below  2ujp  and  uop.  In  this  regime,  we  can  approximate 
7  =  evp/ ptvjJF  where  //  is  the  mobility.  The  ability  to  tune  Uf,  and  consequently  the 
optical  properties,  through  electrostatic  gating  makes  graphene  unique  compared  to 
normal  metals. 

Like  in  noble-metal  plasmonics  (Maier,  2007),  the  free  nature  of  charge  carriers 
described  by  the  Drude  response  gives  rise  to  SP  modes  in  graphene  (Mikhailov  and 
Ziegler,  2007;  Jablan  et  ah,  2009),  which  are  combined  excitations  of  charge-density 
and  electromagnetic  waves  bound  to  the  surface.  At  first  order  in  ksp/kp  the  SP 
dispersion  is  given  by 

e2U]F 

=  2ne  £hksp  ~  2uf  Vf  ksp  (4-1) 

where  £  ~  2.4  is  the  effective  dielectric  constant  and  vp  ~  106  m/s  is  the  Fermi 
velocity  (Wunsch  et  ah,  2006).  This  dispersion  relation  implies  a  remarkable  reduc¬ 
tion  of  the  SP  wavelength  compared  to  the  free  space  wavelength  A0  =  2irc/u>sp,  as 
\sp/ Ao  ~  vf/c  ~  3  •  10-3.  Thus,  the  smallest  possible  mode  volume  of  a  graphene  SP 
resonator,  V  ~  A3p,  can  be  ~  106  times  smaller  than  in  free  space  (Koppens  et  ah, 
2011). 
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4.3  Nonlinear  Plasmonics  in  Graphene 

Because  the  plasmons  are  intra-band  excitations  of  electrons  near  the  Fermi 
surface,  the  nonlinear  conductivity  can  be  calculated  from  the  semiclassical  Maxwell- 
Boltzmann  equations,  as  detailed  in  Appendix  C.l.  In  summary,  the  distribution 
function  f(x,k,t)  for  an  electron  at  in-plane  position  x  and  with  Bloch  momentum 
k  evolves  under  the  Maxwell-Boltzmann  equation  as 


d tf  +  vFk  ■  dxf  +  edx(p  ■  dkf  =  0,  (4.2) 

where  the  electostatic  potential  (f>(x,  z,  t )  satisfies  Poisson’s  equation  V20  =  enS(z)/e oe 
Here  z  is  the  out-of-plane  coordinate  and  n  —  f  dk  f  is  the  2D  electron  density.  For 
weak  excitations  of  the  electron  distribution,  the  term  dkf  in  the  Maxwell-Boltzmann 
equation  can  be  replaced  by  the  equilibrium  value  dkf^°\  yielding  a  linear  equation 
supporting  SPs  with  the  dispersion  given  in  Eq.  (4.1)  and  an  electrostatic  wave  given 
by  E  =  —  V0  oc  8nsm(kx  —  cot). 

For  sufficiently  large  density  perturbations  Sn,  the  nonlinear  interaction  between 
the  non-equilibrium  distribution  dkf  and  potential  must  be  accounted  for.  This  effect 
can  be  interpreted  as  a  backaction  induced  by  the  electrostatic  wave  on  the  electrons 
via  a  ponderomotive  force  Fp  rsj  dxE2  rs_/  kSn 2  sin  2 kx,  which  grows  with  the  amplitude 
of  the  SPs.  This  nonlinear  force  directly  excites  a  second  plasmon  wave  at  wavevector 
2k  and  frequency  2c o,  i.e.  second  harmonic  generation,  and  gives  rise  to  the  second 
order  conductivity  calculated  by  Mikhailov  (2011).  We  show  (see  Appendix  C.l)  that 
this  leads  to  a  nonlinear  shift  at  the  original  wavevector  k  and  frequency  u,  with  an 
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effective  third  order  conductivity  for  the  SPs  given  by 


a(3\ksp,u)  = 


.  3n  Vp  CqE2 
4  ujp  Huj 


(4.3) 


This  result  differs  from  the  nonlinear  conductivity  as  seen  by  free-space  light  normally 
incident  on  a  graphene  sheet,  where  one  finds  that  cd3)  ~  1/cu3  (Mikhailov  and  Ziegler, 
2008).  Remarkably,  as  we  discuss  next,  the  tight  confinement  of  SPs  in  graphene 
implies  that  the  fields  associated  with  even  single  quantized  SPs  are  strong  enough 
that  nonlinear  effects  are  observable. 


4.4  Graphene  Macro- Atom 

Anticipating  the  large  strength  of  nonlinear  interactions  at  the  level  of  single 
SPs  in  nanoscale  graphene  resonators,  we  are  motivated  to  introduce  a  quantum 
description  of  such  a  system.  We  write  the  Hamiltonian  as  H  =  H0  +  Hc ,  where 
Ho  characterizes  solely  the  excitation  spectrum  of  the  graphene  resonator,  and  Hc 
describes  an  external  coupling  to  the  resonator  (such  as  in  Fig.  4.3a),  which  allows 
one  to  probe  the  resonator  properties  or  utilize  the  nonlinearities  for  applications  such 
as  a  single-photon  transistor. 

We  first  consider  the  intrinsic  properties  of  the  resonator  given  by  H0.  Considering 
the  fundamental  SP  mode  of  the  resonator  with  corresponding  annihilation  operator 
aq  and  number  operator  nq  =  a^qaq,  the  effective  Hamiltonian  H0  is  given  by  (Denardo 
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Figure  4.1:  a)  Schematic  of  the  graphene  macro-atom.  A  doped  graphene  disk  confines 
photons  as  plasmons  to  mode  volumes  millions  of  times  smaller  than  free  space.  This 
induces  a  large  dispersive  nonlinearity  so  that  only  a  single  photon  can  resonantly 
excite  the  cavity,  b)  Shows  the  nonlinear  shift  from  Eq.  4.5  for  the  fundamental 
mode  relative  to  the  plasmon  linewidth  with  decreasing  mode  volume  \  q  =  (Asp/A0)3. 
Here  we  took  the  linewidth  as  7  =  evF/phu>F  with  the  Fermi  energy  huiF  =  0.2  eV 
and  a  mobilities  of  p  =  105(104)  cm2 /Vs  corresponding  to  quality  factors  of  roughly 
600(60). 


and  Putterman,  1988;  Gervasoni  and  Arista,  2003) 

Ho  =  {uiq  -  it c/2  +  gq{nq  -  1))  nq.  (4.4) 

See  Appendix  C.2  for  a  detailed  derivation.  This  Hamiltonian  describes  the  quantum 
analog  of  a  cavity  exhibiting  an  intensity-dependent  refractive  index,  where  the  effec¬ 
tive  resonant  frequency  u)q  -f  rjq(nq  —  1)  shifts  depending  on  the  intra-cavity  photon 
number.  Here  we  have  also  included  the  total  cavity  linewidth  k  =  Kex  +  7  into 
the  cavity  description  which  includes  both  the  intrinsic  losses  due  to  impurity  and 
phonon  scattering,  given  by  7  and  radiative  losses  of  the  cavity  into  other  optical  or 
plasmonic-  modes,  given  by  nex.  For  graphene,  the  nonlinear  interaction  strength  is 
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given  by 

7n  u)q  eq3 

T,q  ~  2^’  (4‘5) 

where  ag  =  e2/47reo  hvF  ~  2  and  A  is  the  mode  area  of  the  resonator,  which  can  be 
given  by  A  =  A2p/4  for  a  diffraction-limited  structure.  The  r)q  oc  A^1  scaling  reflects 
that  the  held  intensity  of  a  single  SP  grows  inversely  like  its  confinement. 

At  the  quantum  level,  the  interaction  parameter  2r/g  indicates  the  additional  en¬ 
ergy  cost  to  excite  two  versus  one  photon  in  the  cavity,  as  can  be  seen  in  the  cavity 
excitation  spectrum  (inset  of  Fig.  4.2a).  When  2 rj  k,  the  graphene  sheet  behaves  as 
a  two-level  atom  because  it  can  only  resonantly  absorb  a  single  photon  as  illustrated 
in  Fig.  4.1a;  thus  we  describe  this  as  the  quantum  nonlinear  regime.  The  ratio  2 r]q/ k 
is  then  a  good  measure  of  the  quality  of  the  cavity  as  a  quantum  emitter.  Fig.  4.1b 
shows  2r)q/'y  for  the  fundamental  mode  with  decreasing  mode  volume  (assuming  mo¬ 
bilities  of  105  and  104  cm2/Vs),  where  we  see  that  this  ratio  can  be  as  large  as  100. 
The  parameter  ij/k  oc  Q/A,  where  Q  is  the  quality  factor  of  the  resonator.  Intuitively 
then,  the  nonlinear  cavity  can  exhibit  quantum  effects  if  two  photons  interact  within 
a  small  enough  volume  and  for  long  enough  time. 

The  enabling  mechanism  for  a  two-level  atom  to  be  useful  for  quantum  information 
processing  is  that  it  can  only  emit  single  photons  at  a  time.  This  can  be  characterized 
by  the  second  order  correlation  function  of  the  emitted  light,  which  is  identical  to 
that  of  the  cavity  mode,  g^2\t)  =  (cd(r)at(f  +  r)a(t  +  r))a(r))/(a'l’(r)a(r)).  For  a 
stationary  process,  g^(0)  <  1  indicates  non-classical  “anti-bunching”  and  approaches 
g(2)(0)  =  0  in  the  limit  of  an  ideal  two-level  emitter.  We  consider  the  case  where  the 
resonator  is  driven  by  an  external  laser  from  the  side  and  emission  is  collected  from 
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Figure  4.2:  a)  Shows  g<2>(0)  for  the  graphene  macro-atom  driven  by  a  weak  coherent 
state.  As  the  plasmon  lifetime  increases  g^  (0)  becomes  much  less  than  one,  indicating 
its  transition  to  an  effective  two- level  system  (illustrated  in  inset),  b)  Shows  g{2\t) 
for  hu)sp  =  0.2  eV  and  two  different  mobilities. 


a  different  direction.  In  the  limit  of  weak  driving  we  find  that 


9W(0)  = 


4?72  4-  ac2  ’ 


(4.6) 


thus  establishing  rj  <  k  as  the  regime  where  quantum  properties  become  observable. 
In  Fig.  4.2ab  we  take  nex  =  0  and  we  see  that,  for  the  largest  nonlinearities,  g12^  <  1 
can  be  readily  observed,  even  in  the  presence  of  significant  loss. 


4.5  Efficient  Coupling  and  a  Single-photon  Switch 

In  order  to  exploit  the  large  nonlinearity  of  graphene,  we  need  an  efficient  method 
to  convert  SPs  into  external  optical  modes  on  time  scales  short  compared  to  the 
intrinsic  losses.  Specifically,  one  needs  that  the  total  linewidth  k  =  Kex  -I-  7  contains 
a  large  component  nex  that  goes  into  desirable  external  channels  compared  to  the 
intrinsic  losses  7.  One  approach  is  to  use  the  direct  dipolar  emission  of  the  cavity 
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Figure  4.3:  a)  Integrated  nonlinear  optical  circuit  for  interfacing  the  graphene  macro¬ 
atom  with  photons.  First  the  photons  are  converted  into  bulk  plasmons  via  a  grating, 
then  they  couple  to  the  graphene  macro  atom,  after  which  they  are  converted  back 
into  waveguide  photons,  b)  Shows  a  top  down  view  near  the  plasmon  cavity,  c)  Shows 
the  single  photon  transmission  through  the  system  which  is  less  than  one  due  to  losses 
during  the  grating  coupling,  here  we  took  £  =  ko  and  P  1  so  the  only  losses  are 
in  the  nanoribbons.  The  plasmon  frequency  is  taken  to  be  less  than  the  cutoff  from 
optical  phonons  (0.2  eV)  and  we  assume  the  decay  rate  7  is  dominated  by  impurity 
scattering.  The  three  curves  are  for  a  fixed  plasmon  frequency  with  increasing  Fermi 
energy,  which  increases  the  spatial  propagation  length  of  the  plasmons.  d)  Shows 
bunching  in  reflection  for  two  incident  photons  from  the  left,  with  ftu>sp  =  0.2  eV,  EF  = 
0.23  eV,  P  =  2,  and  mobilities  p  =  104(105)  cnr/V-s  (dashed(solid))  corresponding 
to  a  lifetime  of  0.2(2)  ps  and  a  cavity  quality  factor  of  60(600).  e)  Shows  antibunching 
in  transmission  for  P  —  0.1  with  other  parameters  as  in  (c). 
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into  free  space  radiation.  For  the  square  cavities  described  above,  the  dipole  moment 
is  given  by  p  =  2  e  k2F/ k^p  which  gives  a  decay  rate  into  radiation  of 


kl  p2 

37reo  h 


16  ag  kp 

3  k3sp 


Vq  ojf 


(4.7) 


where  Vo  =  (Asp/Ao)3.  For  cavities  in  the  quantum  nonlinear  regime,  this  is  a  small 
contribution  to  the  total  losses;  thus,  while  it  may  be  a  convenient  method  for  probing 
the  system  a  more  practical  approach  is  needed. 

We  envision  a  two-step  process  illustrated  in  Fig.  4.3ab:  first  a  waveguide  photon 
is  converted  into  a  bulk  plasmon  via  a  dielectric  grating,  then  this  plasmon  can  tunnel 
directly  into  the  nonlinear  cavity.  We  first  consider  the  direct  coupling  between  the 
cavity  and  the  bulk  plasmons.  We  take  the  cavity  to  be  separated  a  distance  d  from  a 
long  nanoribbon  of  width  W.  For  d  W,  Xsp  the  coupling  is  dipolar  and  small,  which 
allows  us  to  calculate  the  decay  of  the  fundamental  cavity  mode  into  the  nanoribbon 
via  Fermi’s  golden  rule  (see  Appendix  C.3) 


32  krF  Wcj 

7t2  PC  b*  H  W6 
n  KspKspa 


(4.8) 


where  kpC  is  the  Fermi  wavevector  in  the  nanoribbon  (r)  and  cavity  (c)  and  k*p  is  the 
wavevector  for  the  nanoribbon  plasmon  that  is  resonant  with  the  cavity  mode.  The 
cavity  can  be  efficiently  controlled  through  the  nanoribbon  by  operating  at  a  distance 
d  such  that  this  decay  is  the  dominant  loss  channel  for  the  cavity. 

Once  the  plasmon  is  in  the  bulk  it  still  remains  to  out-couple  it  to  the  waveguide. 
Due  to  the  large  mismatch  in  wavevectors,  ksp/k0  ~  c/up,  the  bare  coupling  of  the 
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plasmons  to  waveguide  mode  will  be  very  small.  A  simple  and  convenient  solution 
is  to  fabricate  a  dielectric  grating  to  enable  momentum  conservation.  For  parallel 
propagation,  the  grating  wavector  kg  should  be  given  by  kg  =  ksp  —  hp. 

Here  we  consider  the  case  of  a  single-mode  dielectric  slab  waveguide  in  vacuum 
coupled  via  the  grating  to  a  graphene  nanoribbon.  This  geometry  can  be  analyzed 
via  coupled  mode  theory  and  optimized  as  a  function  of  the  slab  thickness  (Snyder 
and  Love,  1983).  Taking  the  grating  profile  to  be  of  the  form  eg(x )  =  Secoskgx  gives 
the  power  conversion  for  weak  losses  between  the  waveguide  and  plasmon  mode  as 
cos2(£a;)  where  £  is  spatial  coupling  between  the  TM  mode  of  the  waveguide  and 
nanoribbon  (see  Appendix  C.4) 

(4.9) 

here  W'  >  W  is  the  width  of  the  waveguide,  y2  =  (32  —  k q  is  the  transverse  wavevector 
of  the  slab  mode,  [3  is  the  longitudinal  wavevector,  and  h  is  the  distance  between  the 
slab  and  the  graphene.  Because  the  factor  in  £  in  front  of  is  order  unity,  the 
plasmon  conversion  for  a  weak  grating  is  limited  to  distances  ~  Ao  3>  A sp.  As  a  result 
the  spatial  decay  rate  of  the  plasmons  must  be  much  larger  than  ko  to  achieve  efficient 
conversion.  For  plasmon  frequencies  below  the  cutoff  from  optical  phonons  (~  0.2 
eV)  the  spatial  decay  rate  is  given  by  'yksp/u>sp  ~  evp  husp/2  p  Ep,  which  strongly 
decreases  with  the  Fermi  energy  Ep.  In  Fig.  4.3c  we  show  the  transmission  of  a  single 
photon  through  the  geometry  displayed  in  Fig.  4.3ab. 

The  device  depicted  in  Fig.  4.3ab  can  be  used  as  a  nonlinear  single-photon  switch. 
To  characterize  this  process,  it  is  first  necessary  to  understand  how  an  input  held 
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through  the  waveguide  is  transformed  upon  interacting  with  the  nonlinear  resonator, 
which  can  be  done  through  an  input-output  formalism.  In  the  case  of  Fig.  4.3ab  of  a 
resonator  equally  coupled  to  two  waveguides,  the  resonator  evolves  under  the  incoming 
fields  of  the  left-  and  right-going  modes  under  the  Hamiltonian  Hc  =  \Jnex{arvn  + 
a\n)a)  +  h.c,,  while  the  output  fields  are  given  by  ar}^t  =  am'1  +  i\[C^,a. 

This  one  dimensional  model  has  been  solved  exactly  for  the  case  of  one  and  two 
resonant  photons  input  from  a  single  direction  in  the  waveguide  (Liao  and  Law, 
2010).  The  response  is  characterized  by  the  effective  Purcell  factor  P  =  nex/^,  which 
measures  the  fraction  of  cavity  emission  into  the  waveguide,  and  the  normalized 
nonlinearity  fj  =  p/n.  The  transmission  t  and  reflection  r  coefficients  for  a  single 
photon  incident  on  resonance  with  the  cavity  are  given  by  t  —  —P/{  1  +  P)  and 
r  =  1/(1  +  P).  The  two  photon  response,  however,  is  modified  by  the  nonlinearity. 
For  example  two  photons  at  frequency  usp  will  be  blocked  from  entering  the  cavity 
due  to  the  nonlinearity.  This  leads  to  antibunching  in  the  transmission  and  bunching 
in  the  reflection  as  shown  in  Figures  4.3de.  The  suppression  in  the  transmission  scales 
as  ’fj2  similarly  to  Eq.  4.6,  while  the  bunching  in  reflection  scales  as  P4  for  fj  P  3>  1 
(Liao  and  Law,  2010).  Fig.  4.3e  shows  that  such  a  device  essentially  realizes  a  single 
photon  transistor  where  one  control  photon  can  block  several  signal  photons  from 
propagating  through  the  cavity  for  a  time  given  by  the  inverse  cavity  lifetime.  In 
addition,  the  signal  photons  do  not  have  to  be  at  the  same  frequency  as  the  control 
photon  so  long  as  they  have  significant  nonlinear  interaction  in  the  cavity. 
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4.6  Conclusions 

Our  analysis  shows  that  graphene  plasmonics  may  provide  a  powerful  platform 
for  the  nonlinear  quantum  optical  control  of  light.  Combined  with  the  scalable  fab¬ 
rication  of  graphene  this  could  allow  the  creation  of  complex  quantum  networks  for 
many  applications  in  quantum  information  and  quantum  simulation,  as  well  as  in 
classical  nonlinear  optics  (Carusotto  and  Ciuti,  2013).  Such  a  system  is  ultimately 
limited  either  by  the  losses  in  graphene  or  the  strength  of  the  nonlinearity.  We  esti¬ 
mate  currently  achievable  quality  factors  for  the  plasmon  cavity  range  from  10  —  103; 
however,  estimates  of  the  ultimate  limit  to  the  graphene  plasmon  lifetime  suggest 
that  quality  factors  greater  than  104  are  possible  (Principi  et  ah,  2013).  To  enhance 
the  nonlinearity  further  hybrid  structures  can  be  envisioned  if  one  could  fabricate  the 
structure  on  top  of  a  strong  nonlinear  substrate. 
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All-Optical  Switch  and  Transistor 
Gated  by  One  Photon 

5.1  Introduction 

In  this  chapter  we  report  on  the  theoretical  analysis  of  an  experimental  real¬ 
ization  of  an  all-optical  transistor  where  one  ‘gate’  photon  controls  a  ‘source’  light 
beam.  Using  a  slowed  a  light  pulse  in  an  atomic  ensemble  contained  inside  an  opti¬ 
cal  resonator,  we  demonstrate  that  one  stored  gate  photon  can  control  the  resonator 
transmission  of  subsequently  applied  source  photons.  In  continuous  operation,  sig¬ 
nal  and  gate  photons  derived  from  different  lasers  become  anti-correlated  with  an 
equal-time  cross-correlation  function  </2)(0)  =  0.89  ±  0.01. 

Photons  are  excellent  carriers  of  quantum  information,  but  it  is  difficult  to  in¬ 
duce  the  strong  interactions  between  individual  photons  that  are  required  for,  e.g., 
all-optical  quantum  information  processing.  Nevertheless,  advances  toward  such  in- 
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teractions  have  been  made  in  cavity  quantum  electrodynamics  (QED)  systems  with 
atoms  (Birnbaum  et  al.,  2005;  Brennecke  et  ah,  2007;  Colombe  et  al.,  2007;  Kubanek 
et  al.,  2008;  Tanji-Suzuki  et  al.,  2011;  Brooks  et  al.,  2012)  or  artificial  atoms  (Michler 
et  ah,  2000;  Press  et  ah,  2007;  Fushman  et  ah,  2008;  Volz  et  ah,  2012;  Bose  et  ah,  2012), 
and  in  a  cavity-free  system  using  atomic  Rydberg  states  (Dudin  and  Kuzmich,  2012; 
Peyroncl  et  ah,  2012)  or  dye  molecules  (Hwang  et  ah,  2009).  All-optical  switching  of 
one  beam  by  another(Bajcsy  et  ah,  2009)  and  cross-phase  modulation(Lo  et  ah,  2011) 
have  been  demonstrated  at  the  level  of  a  few  hundred  photons  by  means  of  electromag- 
netically  induced  transparency  (EIT)(Flcischhauer  et  ah,  2005).  At  the  few-photon 
level,  nonclassical  light  has  been  generated  (Michler  et  ah,  2000;  Birnbaum  et  ah, 
2005;  Press  et  ah,  2007;  Kubanek  et  ah,  2008;  Dayan  et  ah,  2008;  Fushman  et  ah, 
2008;  Bose  et  ah,  2012;  Dudin  and  Kuzmich,  2012;  Peyroncl  et  ah,  2012;  Brooks  et  ah, 
2012),  and  optical  nonlinearities  of  16°  in  phase  shift (Turchette  et  ah,  1995)  and  up 
to  ~20%  in  two-photon  attenuation  (Fushman  et  ah,  2008;  Tanji-Suzuki  et  ah,  2011; 
Volz  et  ah,  2012)  have  been  observed  in  cavity  QED  systems.  While  switching  of  the 
cavity  transmission  by  a  single  atom  has  also  been  achieved  (Thompson  et  ah,  1992), 
the  realization  of  an  optical  transistor  exhibiting  gain  with  gate  signals  at  the  few- 
or  one-photon  level  (Chang  et  ah,  2007)  remains  a  challenge. 

Here  we  theoretically  analyze  a  cavity  QED  version  (Grangier  et  ah,  1998;  Imamoglu 
et  ah,  1997)  of  a  single-photon  switch  based  on  EIT  in  a  four-level  system  (Fleis- 
chhauer  et  ah,  2005).  Combining  this  technique  with  slow  light  allows  one  to  im¬ 
plement  an  all-optical  transistor  where  one  gate  photon  can  switch  multiple  signal 
photons.  The  device  performance  is  quantified  by  measuring  interaction-induced 
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photon-photon  anticorrelations  between  two  distinct  modes  driven  by  independent 
lasers. 

The  system  (Tanji-Suzuki  et  al.,  2011)  consists  of  an  ensemble  of  laser-cooled  ce¬ 
sium  atoms  optically  trapped  inside  a  high-finesse  optical  cavity  (Fig.  5.1  A)  operating 
in  the  strong-coupling  regime  (Birnbaum  et  ah,  2005;  Brennecke  et  ah,  2007;  Colombe 
et  ah,  2007;  Kubanek  et  ah,  2008;  Tanji-Suzuki  et  ah,  2011;  Brooks  et  ah,  2012)  of 
cavity  QED.  Each  atom  has  a  four-state  IV- type  level  structure  | g)  «->•  |  d)  O  |s)  |e) 

with  two  stable  ground  states  j g),  s),  and  two  electronic  excited  states  | d),  \e) 
(Fig.  5. IB).  For  atoms  prepared  in  state  | g),  this  atomic  structure  mediates  an  ef¬ 
fective  interaction  between  free-space  photons  (photons  resonant  with  the  | g)  — >  \d) 
transition  serving  as  gate  photons)  and  cavity  photons  (photons  resonant  with  the 
|s)  — >  \e)  transition  serving  as  the  source) (Schmidt  and  Imamoglu,  1996;  Imamoglu 
et  ah,  1997;  Harris  and  Yamamoto,  1998).  These  two  transitions  are  connected  via  a 
control  laser  that  addresses  the  |  d)  — *  s)  transition  and  induces  transparency  (E1T) 
for  the  gate  photons. 

Without  the  signal  beam,  the  gate  photons  are  transmitted  through  the  ensem¬ 
ble,  traveling  in  the  medium  as  slow-light  polaritons,  a  superposition  of  a  photon 
and  a  collective  atomic  excitation  to  the  state  s) .  In  the  absence  of  gate  photons, 
the  state  |s)  is  unpopulated  so  resonant  signal  photons  are  transmitted  through  the 
cavity.  Thus  photons  arriving  individually  in  either  the  gate  or  the  signal  mode  are 
transmitted  through  the  system.  However,  when  photons  are  simultaneously  present 
in  the  two  modes,  they  affect  each  others  propagation.  A  signal  photon  inside  the 
cavity  introduces  a  decoherence  path  for  the  state  s)  via  coupling  to  the  unstable 
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excited  state  \e)  and  reduces  the  EIT  transmission  for  the  gate  photon.  Conversely, 
when  a  gate  photon  travels  in  the  atomic  medium  as  a  slow-light  polariton,  the  po- 
laritons  atomic  component  in  state  .s)  reduces  the  cavity  transmission  by  introducing 
additional  cavity  loss  through  photon  scattering  on  the  |s)  — >■  |e)  transition.  If  the 
coupling  between  the  cavity  and  a  single  atom  in  state  |s)  is  sufficiently  strong,  a  sin¬ 
gle  photon  in  the  gate  or  the  signal  mode  will  each  block  the  other  mode  (Thompson 
et  al.,  1992). 

The  strength  of  this  effective  photon-photon  interaction  essentially  depends  on 
two  parameters:  the  free-space  resonant  optical  depth  J\f  on  the  gate  transition  that 
measures  the  collective  coupling  of  the  atomic  ensemble  to  the  gate  photon,  and 
the  single-atom  cavity  cooperativity  r]  that  sets  the  interaction  strength  between  one 
atom  in  |s)  and  one  signal  photon  (Birnbaum  et  ah,  2005;  Kubanek  et  ah,  2008;  Tanji- 
Suzuki  et  ah,  2011).  The  optical  depth  A f  sets  an  upper  limit  of  1  —  e~^  for  the  atomic 
component  |s)  of  the  slow  light  polariton  (Fleischhauer  et  ah,  2005)  (the  component 
that  can  interact  with  the  cavity  mode).  The  reduction  of  the  cavity  transmission 
by  a  single  intracavity  atom  in  state  |s),  on  the  other  hand,  is  given  (Tanji-Suzuki 
et  ah,  2011)  by  the  factor  (1  +  r/)-2.  For  deterministic  two-mode  photon-photon 
interactions,  the  gate  photon  must  be  converted  with  reasonable  efficiency  into  an 
atomic  population  in  state  |s)  which  then  needs  to  block  the  cavity.  Hence  we  require 
both  strong  collective  coupling  on  the  gate  transition  (A/"  A>  1)  and  strong  singlc- 
atom-cavity  coupling  (Imamoglu  et  ah,  1997;  Grangier  et  ah,  1998)  on  the  signal 
transition  (?/  3>  1). 
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Figure  5.1:  All-optical  single-photon  transistor.  Setup  (A)  and  atomic  level  scheme 
(B).  An  ensemble  of  laser  cooled  atoms  is  trapped  inside  an  optical  resonator  operating 
in  the  single-  atom  strong-coupling  regime.  The  atoms  are  prepared  in  state  | g) 
by  optical  pumping,  and  the  coupling  beam  on  the  |s)  — >  \d)  transition  induces 
transparency  (EIT)  for  gate  photons  on  the  \g)  — »  \d)  transition.  The  gate  photons 
travel  slowly  through  the  medium  as  quasi-  particles  (dark-state  polaritons)  with  an 
atomic  spin-excitation  component  in  the  state  |s)  that  interacts  with  signal  photons 
on  the  |s)  -*  |e)  transition.  The  interaction  results  in  photon-photon  anticorrelations 
that  are  measured  with  photon  counters  Dg  and  Ds.  The  atomic  states  of  133Cs 
used  in  this  experiment  are  \y)  =  1 65^/2-.  F  =  3,  mF  =  3)  ,  |d)  =  |6P3/2, 4, 4)  ,  |s)  = 
|6S'i/2i  4.4) ,  \e)  —  |6P}/2, 5, 5),  where  F  and  tuf  denote  the  hyperfine  and  magnetic 
sublevels. 
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5.2  Theoretical  Model 

Including  the  decay,  the  effective  Hamiltonian  for  this  system  can  be  written  as 
Wf,/h  =  £  c \k\  a\ak  +  (uc  —  in/2)Cb  +  i£(tf  —  b) 

k 

+  (ugd  -  5]  \d)i  (d\  +  (use  -  1 e)i  (el  (5-1) 

i  i 

+  {Qc  elUdst  \s)x  (d\  +  g9  a\x)  \g)x  (d\  +  gs  |s)i  (e|  +  h.c .) 

i 

Here,  c  is  the  speed  of  light,  k  is  the  wavenumber  of  the  gate  held,  ojc  is  the  cavity 
frequency  and  n  is  the  decay  rate  of  the  cavity.  The  electric  held  operators  for  the 
two  helds  can  be  written  as  £g(x)  =  \J^^a(x)  and  £s  =  sj where  a{x )  = 
N~ 1/2  elkxak  and  b  are  bosonic  annihilation  operators,  cko  is  the  center  frequency 
of  the  gate  held,  and  V  is  the  quantization  volume.  Additionally,  £  is  the  amplitude 
of  the  cavity  input  held,  c is  the  atomic  transition  energy  between  states  g  and  z/, 
If  is  the  classical  Rabi  held  for  the  coupling  held,  T  is  the  linewidth  of  the  excited 
states  | d)  and  |e),  7  is  decoherence  rate  of  two  stable  ground  states  | g)  and  s),  and 
gg,  gs  are  the  bare  couplings  of  the  atomic  transition  to  the  two  helds.  We  take  the 
gate  and  signal  helds  to  be  resonant  with  the  atoms  so  that  cko  =  w gd  and  uic  =  00 se. 

The  use  of  this  effective  Hamiltonian  is  sufficient  to  describe  the  steady  state  for 
the  case  of  weak  input  helds  gg  (cba)  <C  02/T  and  gs  (b^b)  k.  In  this  limit  we  can 
take  the  approach  of  Carmichael  et  al.  (1991)  to  calculate  the  two  time  correlation 
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function  between  the  two  fields 


t  +  r)a(x ,  t  +  r)b(t)) 
(at (x,  t)a(x ,  t))(&t ( t)b(t )) 


(5.2) 


In  this  limit  we  can  write  the  density  matrix  as  a  product  state  p  =  |x(t))  (x(r)| 
and  we  truncate  the  available  states  in  the  system  at  the  level  of  two  excitations  from 
the  state  with  zero  photons  and  all  atoms  in  | g),  which  we  refer  to  as  \g,  0,  0).  The 
one  excitation  states  are  \g,  lx,0)  =  a^(x)  \g,  0,  0),  |g,0,l)  =  6^  |gr,  0,  0) ,  |dx,0,  0)  = 
<jjg  \g,  0,  0),  and  |sx,0,  0)  =  cr*g  \g,  0,  0),  where  =  \g)x(v\.  The  two  excitation 
states  that  are  relevant  for  gg2J  are  |g,lx,l)  =  a^(x)  |g,0, 1),  |dx,0, 1)  =  b'  \dx,0,  0), 
|s*,  0, 1)  =  b]  |sx,  0,  0),  and  \ex,  0,  0)  =  a*s  \sx,  0,  0). 

We  then  expand  \x(t))  in  these  states  and  find  the  evolution  according  to  = 
Heff  |x)  applying  the  boundary  condition  that  the  free  space  input  field  is  a  weak 
coherent  state.  The  only  terms  in  Heff  which  create  excitations  are  the  driving  fields, 
which  are  perturbative  implying  that  the  amplitude  of  the  one  excitation  states  are 
proportional  to  £  and  the  two  excitation  amplitudes  are  proportional  to  £2. 

To  calculate  gf}(r)  we  take  the  picture  where  the  detection  corresponds  to  a 
quantum  jump  from  the  steady  state  |xss)  into  the  state  a(x,t )  |xss)  for  r  <  0  and 
b(t)  |xss)  for  r  >  0  (Carmichael  et  ah,  1991).  To  find  ggs(r)  we  can  then  simply 
evolve  the  operator  ns(t)  or  ng(x,t )  for  a  time  r  under  Heff  starting  from  the  jump 
state. 

To  find  the  steady  state  we  expand  |x(i))  in  the  zero,  one  and  two  excitation 


79 


98 


Chapter  5:  All- Optical  Switch  and  Transistor  Gated  by  One  Photon 


states 

1x0*5  t))  =  I/,  0, 0)  +  Al(x)  I/,  lx,  0)  +  A\(x)  I gx,  0, 0)  +  A\{x)  \dx,  0, 0)  +  Al \f,  0, 1) 
+  A\(x)  I/,  lx,  1)  +  A%(x)  | gx,  0, 1)  +  Al(x)  \dx,  0, 1)  +  A\(x)  \ex,  0, 0) 

(5.3) 

where  we  neglect  the  one  and  two  excitations  in  the  normalization  because  they  are 


perturbative.  The  equations  of  motion  are  for  the  Aj  are  found  from  i^-  =  Heff  y). 

{d t  +  cdx)Al(x)  =  - igfdy/N  A\(x),  (5.4) 

dtA\(x)  =  — T /2A^{x)  —  iggy/N  A^(x)  —  iOl/2A\(x),  (5.5) 

dtA\(x)  =  —7/2  A\{x)  —  iOl/2A\(x),  (5.6) 

dtA\  =  -k/2A\  +  £,  (5.7) 

(dt  +  cdx)A\(x )  =  —k/2A\(x)  +  SA](x)  -  iggVNAl(x),  (5.8) 

dtA\(x)  =  —  (T  +  k)/2AI(x)  —  iggVN A\(x)  —  ifl/2Al(x)  +  SAl(x),  (5.9) 

dtAl(x)  =  ~(k  +  ry)/2Al(x)  —  igsAA(x )  —  iVl/2A\{x)  +  £A\(x),  (5.10) 

dtAl(x)  =  -T/2A\(x)  -igsA\(x)  (5.11) 
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(2\ 

These  eight  equations  are  the  only  ones  relevant  for  ggS  (t),  they  give  the  steady  state 


2  glN  x 


Al(x)  =  a  exp  f  — 


r  +  o2/7c 


=  a  exp  l  — 


A f 


x 


2(1  +  fi2/7r)  l 


Al(x)  = 


fi2  +  7r 

i2gg\/N 

~r  +  n2/7 


Al{x), 


A\  =  ^~, 

3  k/2’ 


V 


Al(x)  1 

Al(x)Al  1  +  rj^  1  +  r) 


exp 


~U]+o(k/t) 


(5.12) 

(5.13) 

(5.14) 

(5.15) 

(5.16) 


where  a  is  the  amplitude  of  the  input  coherent  state,  N  is  the  number  of  atoms, 
J\f  =  Ag2gNL/cT  is  the  optical  depth,  L  is  the  length  of  the  medium,  r)  =  4g2 / kP  is 
the  cooperativity,  and  we  have  defined 


c 


/  +  ny«r  +  7 AA 


(5.17) 


a  correction  factor  that  arises  from  imperfect  E1T,  which  reduces  the  effective  switch¬ 
ing  by  decreasing  the  atomic-excitation  component  of  the  polariton. 

When  r  <  0  the  free  space  photon  is  detected  first  leading  to  a  quantum  jump 
into  the  state 


I  Xj) 


_ ajV)  |y88) _ 

V(Xss\al{L,T)a(L,T)  |x«) 


I/,  0,0)  + 


I/,  0,1) 


(5.18) 


Now 
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(2)/_\  =  {Xj{t)\bHT)b{T)\\j(t))  [ij(Ql2 
%A  '  (x..|M(r)6(r)|x„>  |A>|2 

=  l-(l-  ' 

[  V  >  l  +  <7  J 


(5.19) 

(5.20) 


where  k<  =  n.  For  r  >  0  the  procedure  is  the  same,  except  we  have  to  evolve 
Eqs.  4-6  starting  from  the  initial  conditions  Aq(x,t)  =  A\(x),  A\(x,t)  =  A%(x),  and 
A\(x,t )  =  A\{x).  This  corresponds  to  the  state  |xj)  oc  6|xss)-  The  result  can  be 
expressed  in  the  same  form  as  Eq.  5.19  with  k<  replaced  by  k>  =  Q2/T  +  7 (vg/L)  in 
the  limit  of  small(large)  optical  depth. 


/o\ 

Figure  5.2:  Time  ordered  cross-correlation  function  ggs  (t)  of  the  gate  and  signal  field 
in  steady  state  in  the  limit  of  large  A f.  Parameters  are  such  that  the  cooperativity  is 
7]  =  gg/nT  =  4.5,  the  group  velocity  for  the  gate  field  is  vg/c  =  0*/ g\N  =  10-4,  and 
k  =  10_3c/L. 
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5.3  Experimental  Results 

In  the  limit  of  large  optical  depth  M/C  3>  In  rj  every  incoming  gate  photon  is 
converted  into  a  slow-light  polariton  (Fleischhauer  et  ah,  2005)  with  a  near-unity 
atomic  excitation  component  in  state  |s).  In  this  case,  the  cavity  transmission  is 
modified  by  the  single  intracavity  atom  (Tanji-Suzuki  et  ah,  2011)  and  is  given  by 
(1  +  7/)~2,  which  is  also  the  minimum  value  of  the  cross-correlation  function  g^(0) 
in  this  limit.  This  is  the  limit  shown  theoretically  in  Fig.  5.2  In  the  opposite  limit  of 
large  cooperativity  g  S>  1  and  moderate  optical  depth,  M/C  In  //,  the  signal  photon 
completely  destroys  EIT,  and  gg/  (0)  «  e-Ar  is  simply  the  probability  for  the  gate 
photon  to  pass  through  the  absorbing  medium  in  the  absence  of  EIT.  Interestingly, 
the  correlation  function  g///  (t)  is  asymmetric  in  the  time  separation  t  between  the 
photons  (Hennessy  et  ah,  2007).  This  can  be  understood  as  follows:  the  detection 
of  a  signal  photon  at  time  t  =  0  implies  that  the  EIT  transmission  must  have  been 
reduced  for  times  t  <  0  on  a  time  scale  on  the  order  of  the  cavity  lifetime  k-1,  and 
will  approach  its  uncorrelated  steady-state  value  ggs  =  1  for  times  t  >  0  with  a  time 
constant  determined  by  the  polariton  lifetime,  which  depends  on  the  EIT  linewidth 
02/r  in  the  limit  of  small  optical  depth.  (An  analogous  argument  can  be  made  if  one 
assumes  the  gate  photon  to  be  detected  at  t  =  0.) 

The  measured  cross-correlation  function  ggs (t)  displayed  in  Fig.  5.3  shows  that 
photons  in  the  two  modes  are  uncorrelated  for  large  time  separation  t,  but  display 
a  marked  anticorrelation  dip  near  t  =  0:  when  the  two  photons  derived  from  inde¬ 
pendent  lasers  arrive  near-simultaneously,  they  reduce  each  others  transmission.  The 
data  in  Fig.  5.3  are  well  described  by  Eq.  5.19  using  a  three-parameter  fit  with  the 
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-4  -2  0  2  4 

T  =  ig-ts(tys) 


Figure  5.3:  Mutual  photon-photon  switching  in  continuous  operation.  The  second- 

/  o\ 

order  cross  correlation  function  ggS  (t)  is  displayed  versus  time  separation  r  =  tg  — 
ts  between  photons  in  the  gate  and  signal  modes  for  different  coupling  beam  Rabi 
frequencies:  (A)  TI/2tt  =  0.5  MHz,  (B)  fl/2ii  =  0.9  MHz,  and  (C)  H/27T  =  1.5  MHz. 
Fits  of  the  data  to  the  model  (see  text)  yield  (A)  gg2J( 0)  =  0.91  ±  0.01,  k<  =  (1.6  ± 
0.2)  fis- 1,«>  =  (1.4T0.2)  /is-1,  (B)  ^(0)  =  0.89±0.01,k<  =  (1.5T0.1)  /rs-1,^  = 
(2.4  ±  0.3)  /is-1,  and  (C)  =  0.90  ±  0.01,  «<  =  (1.6  ±  0.1)  /is-1,K>  =  (3.9  ± 

0.3)  /ts-1.  (D)  Fitted  rate  constants  n><  versus  EIT  linewidth  Vt2 /T.  The  positive¬ 
time  rate  constant  fits  to  =  aVt2 /T  +  b  with  slope  a  =  1.0  ±  0.1  and  y-axis 
intercept  b  =  27t(190  ±  30)  kHz,  which  agrees  with  the  expected  values  a  =  1  and 
b  =  7  =  27r(179  ±  10)  kHz.  The  negative-time  rate  constant  tt<  =  1.6  /ts-1  is 
independent  of  the  coupling  laser  intensity  and  is  larger  than  the  cavity  linewidth 
k  =  0.89  /is-1  due  to  incomplete  optical  pumping  leading  to  cavity  absorption.  The 
measurements  were  performed  at  photon  numbers  (ns)  ~  0.1  and  (ng)  ~  0.2  when 
integrated  over  a  time  windows  1/k<  and  1/k>,  respectively. 
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/  o\ 

zero-time  value  gV  (0)  and  the  two  decay  rate  constants  ft/7.  The  time  constants 
obtained  from  the  fit  confirm  the  asymmetric  shape  of  the  cross-correlation  function 
(Fig.  5. 3D):  ft<  =  (1.6  ±0.1)  /xs_1  is  independent  of  control  Rabi  frequency  and  80% 
larger  than  the  cavity  linewidth  ft,  presumably  due  to  the  occasional  presence  of  ab¬ 
sorbing  atoms  in  the  state  |s)  due  to  imperfect  optical  pumping.  (Note  that  even  only 
one  atom  out  of  N  =  2  ■  104  is  sufficient  to  substantially  increase  the  cavity  linewidth 
by  a  factor  I  +  77.)  On  the  other  hand,  the  positive-time  constant  ft>  is  linearly  depen¬ 
dent  on  coupling  beam  intensity  and  agrees  with  the  prediction  D2/T  ±7,  where  fl2/T 
and  7/(271)  =  (179  ±  10)  kHz  have  been  independently  determined  from  separately 
measured  EIT  spectra.  The  fitted  ggf(0)  is  between  0.89  ±0.01  and  0.91  ±0.01  for  the 
three  values  of  the  coupling  Rabi  frequency.  This  agrees  well  with  the  prediction  from 
Eq.  5.19  with  values  between  0.87  ±  0.02  and  0.93  ±  0.02,  using  the  independently 
measured  optical  depth,  fitted  ft><,  and  reduced  cooperativity  7/  =  tik/k<  due  to 
imperfect  optical  pumping.  The  measured  zero-time  correlation  gg2s\ 0)  corresponds 
to  a  mutual  photon-photon  switching  efficiency  of  1  —  gKgs  (0)  =  11%  of  one  photon 
by  the  other  in  continuous  operation. 

5.4  Conclusions 

This  system  constitutes  a  testbed  in  which  we  have  explored  the  physical  principles 
relevant  to  an  all-optical  transistor  based  on  cavity  QED  with  an  atomic  ensemble. 
Before  it  can  be  used  as  a  practical  device,  it  will  be  necessary  to  improve  the  input 
and  output  coupling  efficiencies  for  the  gate  and  source  photons.  The  combined 
storage  and  retrieval  efficiency  of  3%  for  the  gate  photon  is  limited  primarily  by  the 
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optical  density.  The  latter  could  be  improved  by  using  a  deeper  trap,  in  combination 
with  further  cooling  of  the  atomic  ensemble,  which  would  also  increase  the  gate  photon 
storage  time  that  is  currently  limited  by  Doppler  broadening.  The  cavity  outcoupling 
efficiency  for  the  source  photons  of  0.66  could  be  improved  to  0.97  by  using  state- 
of-the-art  mirrors  (Birnbaum  et  ah,  2005;  Brennecke  et  al.,  2007;  Kubanek  et  al., 
2008). 

The  present  work  opens  up  new  perspectives  for  all-optical  information  processing 
with  strong  deterministic  interactions  between  initially  uncorrelated,  distinguishable 
photons.  The  correlations  between  one  gate  and  multiple  source  photons  produced 
by  the  effective  photon-photon  interaction  can  be  used  to  create  two-mode  entangled 
states  of  many  photons.  Finally,  cavities  with  larger  cooperativity  (Birnbaum  et  ah, 
2005;  Brennecke  et  ah,  2007;  Colombe  et  ah,  2007;  Kubanek  et  ah,  2008),  may  enable 
high-fidelity  deterministic  photonic  quantum  gates. 
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Chapter  6 

Few  Body  Physics  in  Strongly 
Interacting  Rydberg-Polariton 
Gases 

6.1  Introduction 

Electromagnetically  Induced  Transparence  (EIT)  can  be  used  to  convert  photons 
coherently  into  atomic  excitations  and  back  (Fleischhauer  et  ah,  2005).  However, 
EIT  by  itself  is  linear  in  the  photon  field,  and  as  such  cannot  be  used  to  induce 
interactions  between  individual  photons.  In  the  optical  transistor  described  in  the 
previous  chapter,  the  gate  photon  is  converted  into  an  atomic  excitation  by  means  of 
EIT,  and  then  the  strong  interaction  between  the  one  excited  atom  and  the  source 
photons  is  accomplished  via  cavity  QED,  resulting  in  strong  coupling.  In  an  alter¬ 
native  free-space  approach  one  can  directly  realize  strong  photon-photon  interactions 
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via  EIT  involving  atomic  Rydberg  levels  with  strong  mutual  atom-atom  interactions. 
The  basic  idea  is  that  while  a  photon  is  traveling  through  the  medium  as  a  Rydberg 
polariton  with  substantial  population  amplitude  in  the  Rydberg  level,  within  a  cer¬ 
tain  characteristic  distance  range  of  the  Erst  photon  (the  so-called  blockade  radius) 
the  second  photon  cannot  experience  EIT  because  the  Rydberg  level  is  shifted  due  to 
the  atomic  Rydberg-Rydberg  interaction  (Lukin  et  ah,  2001;  Pritchard  et  ah,  2010; 
Gorshkov  et  ah,  2011).  For  a  sufficiently  dense  atomic  sample,  such  that  a  photon 
can  be  absorbed  on  a  distance  scale  comparable  to  the  blockade  radius,  this  causes 
optical  nonlinearities  at  the  level  of  individual  photons. 

The  basic  mechanism  underlying  the  interaction  is  outlined  in  Fig.  6.1a.  The 
probe  light  £p  couples  the  ground  state  | g)  to  the  Rydberg  state  |r)  via  an  unstable 
intermediate  state  \e)  of  linewidth  T/ 2  by  means  of  a  control  held  that  is  detuned 
below  the  resonance  frequency  of  the  upper  transition  |e)  — >  |r)  by  A  (Fig.  6.1a.  Un¬ 
der  these  conditions,  EIT  is  established  when  the  probe  detuning  matches  that  of  the 
control  held.  However,  the  Rydberg  medium  is  extremely  nonlinear  and  the  medium 
quickly  saturates  due  to  the  Rydberg  blockade  (Gorshkov  et  al.,  2011;  Pritchard  et  ah, 
2010).  This  results  in  a  two  photon  spectrum  close  to  the  bare  two  level  response, 
such  that  when  A>T  the  nonlinearity  is  purely  dispersive. 

Previous  work  has  demonstrated  nonlinear  phase  shifts  at  the  level  of  two  photons 
and  showed  that  these  phase  shifts  are  associated  with  the  formation  a  two  photon 
bound  state  in  the  steady  state  response  (Peyroncl  et  al.,  2013).  Fig.  6.1b  shows  an 
example  of  such  a  state,  as  evidenced  by  the  bunching  in  the  two  time  correlation 
function  g(-2\t2  —  U)  that  appears  when  the  time  separation  between  the  two  photons 
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Figure  6.1:  a)  Schematic  of  the  Rydberg  mechanism  leading  to  large  dispersive  non- 
linearities.  rs  is  the  blockade  radius,  i\  is  the  EIT  control  field,  £p  is  the  probe  field, 
T  is  the  linewidth  of  the  intermediate  state  |e)  and  Vint  is  the  Rydberg  interaction 
potential,  b)  Numerical  simulation  of  g^2{t)  in  steady  state  for  weak  probe  fields.  We 
took  A/r  =  2  a  Rydberg  blockade  radius  of  10  pm,  a  medium  of  length  100  pm, 
flc  =  10  Mhz  and  an  optical  depth  of  15,  corresponding  to  a  slow  fight  group  velocity 
of  vg  =  600  m/s. 


goes  to  zero.  This  figure  was  obtained  from  numerical  simulations  of  the  steady 
state  two  photon  problem  similar  to  what  is  shown  in  (Peyronel  et  al.,  2013).  In  what 
follows  we  look  more  closely  at  the  time  dependent  dynamics  of  the  two  photon  states, 
including  their  formation  and  the  change  in  their  group  velocity  as  the  interaction 
strength  is  increased.  We  then  look  at  the  formation  of  three  photon  bound  states  in 
steady  state  and  calculate  the  three  time  correlation  function  g^(ti,  t2,  t3). 
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6.2  Numerical  Approach  for  Atom-Photon  Inter¬ 
actions 

Atoms  are  nonlinear  optical  elements  that  can  be  used  to  make  single  photon 
cavities,  photon  transistors,  etc.  Simulation  of  even  single  atom-single  photon  inter¬ 
action  problem  is  non-trivial,  particularly  if  photon  is  not  a  ”  single-mode”  (e.g.  lives 
in  a  waveguide).  Our  approach  is  to  directly  simulate  the  unitary  evolution  of  the 
atom(s)-photon(s)  system  wavefunction, 

|T(f))  =  eim|T(0))  (6.1) 

where  the  Hamiltonian  is 


H  =  H0  +  H' 

(6.2) 

tfo  =  £  ekakak  +  +  -*•)’ 

k  a 

(6.3) 

H'  =  S^J9a(alaa-  +  arna+). 

(6.4) 

a 


The  complication  is  that  photon  parts  of  Hq  and  H'  are  diagonal  in  different 
spaces,  in  momentum  and  real,  respectively.  To  deal  with  this  complication  we  can 
use  the  Trotter  decomposition: 
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| if}{t))  =  re~i^dt'H\^))  =  Uie-iHAt\^(0))  «  nie-</f#Ate-<HoAt|V’(0))  (6.5) 
=  n,:  ^\r)(r\e~lH'At\r)(r\k)(k\e~lHoAt\k)(k\^  |^(0)).  (6.6) 

In  other  words,  the  state  at  a  given  time  step  is  expressed  in  terms  of  the  state 
on  the  previous  step  as 

Mt,))  =  e~iH'At\r)(r\k)(k\e-iHoAt\k)(k\ip(ti-i)).  (6.7) 

By  (Fourier)  transforming  between  momentum  and  real  space  bases,  each  step  of 
evolution  can  be  processed  very  efficiently.  Here  we  apply  these  ideas  to  the  simulation 
of  photons  and  strongly  interacting  atoms  in  Rydberg-Polariton  systems. 

6.3  Results 

Following  Gorshkov  et  al.  (2011)  we  let  £j,(z),  ^(z)  and  V\z)  be  the  slowly 
varying  operators  for  the  creation  of  a  photon,  a  Rydberg  state  |r)  and  an  intermediate 
state  |e),  respectively.  They  satisfy  the  commutation  relations  [£(z),  £{z')\  =  5(z—z'), 
[<S(z),  «S(^,)j  =  6(z  —  z')  and  [P(z),V(z')\  =  S(z  —  z').  The  equations  of  motion  are 


dt£p  =  -cdz£  +  igv/2V 

(6.8) 

dtV  =  -T/2V  +  igp/2£  +  iklc/2S 

(6.9) 

dtS  =  inc/2V  -  i  J  dz'V(z  -  z')S\z')S(z’)S(z) 

(6.10) 
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where  gp  is  the  collective  atom  photon  coupling  which  we  define  by  g2p  =  Tc/£a  where 
Ia  =  L/OD  is  the  absorption  length  and  OD  is  the  optical  depth.  Viz  —  z ')  = 
Cq/(z  —  z'Y  is  the  Rydberg  interaction  which  is  characterized  by  the  blockade  radius 
7'b  =  (2C6r /Ctc)2  defining  the  boundary  where  V (r)  is  greater  than  the  EIT  linewidth 

ayv. 

By  projecting  Eqs.  6.8-10  onto  the  two  excitation  manifold  we  obtain  a  closed 
set  of  equations  for  the  two  photon  dynamics,  which  we  solve  numerically  using 
the  techniques  described  in  section  6.2  As  an  example  Fig.  6.2a  shows  a  pulse  with 
detuning  A  =  2r  and  no  two  photon  detuning  after  traveling  through  the  medium. 
We  clearly  see  the  bunching,  indicating  the  presence  of  the  bound  state  analogously 
to  Fig.  6.1b. 

Since  the  Rydberg  interaction  is  a  short  range  interaction  (note:  a  1  /rn  potentials 
is  considered  short  ranged  if  n^d  the  dimension)  it  is  reasonable  to  ask  whether  it  can 
be  approximated  by  a  delta  potential.  Therefore,  as  an  ansatz  for  the  system  we  use 
a  modified  Nonlinear  Schrodinger  Equation  (NLSE)  governed  by  the  Hamiltonian 


H 


+ 


1  d 2 

2m  dx2 


U 


6(Xi  ~ 


Xi 


hj 


(6.11) 


where  vg  =  c/(l  P  g2/Vt2)  is  the  EIT  group  velocity  and  m  =  —  ^  j-  ^  is  the 
effective  photon  mass  arising  from  the  finite  bandwidth  of  EIT  (Fleischhauer  et  ah, 
2005).  The  interaction  paramter  U  is  given  by 


=  c°  r2  0D b 
8  m  A2  I2a  8 m  A2  la 


(6.12) 
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where  0DB  =  rs/H-a  and  Cq  =  0.92  is  a  numerical  factor.  In  order  for  the  NLSE  to 
be  self-consistent  it  is  necessary  that  the  extent  of  the  bound  states  is  much  larger 
than  the  blockade  radius.  Note  that  U  is  negative,  which  normally  corresponds  to 
repulsion,  but  because  m  <  0  it  acts  as  an  attractive  potential. 


6.3.1  Two  Photon  Solitons 

Solving  for  the  two  photon  bound  states  with  Eq.  6.11  gives  the  dispersion  of  the 
bound  states  as  (Ben-Aryeh,  1999) 


EB(k)  =  vgk+^~ 
2  m 


-mU2 

4 


(6.13) 


To  calculate  the  expected  change  in  group  velocity  in  this  model  we  note  that  since 
our  input  state  is  at  zero  two  photon  detuning  it  excites  the  bound  state  at  EB(k)  =  0 
which  is  shifted  away  from  k  —  0.  This  gives  rise  to  additional  phase  accumulation 
as  discussed  by  Peyronel  et  al.  (2013),  but  also  a  change  in  the  group  velocity 


2  A2 


OD%-  1 


(6.14) 


Since  the  approximation  of  a  delta  potential  is  uncontrolled  we  also  calculated 
the  change  in  vg  using  full  time  dependent  simulations.  In  Fig.  6.2  these  results  are 
compared  to  Eq.  6.14  where  we  see  that  the  agreement  is  good  for  small  ODb,  but 
they  diverge  at  larger  ODB  suggesting  other  effects  are  becoming  important  and  the 
NLSE  is  not  sufficient  to  describe  the  dynamics.  Finally,  we  remark  that  the  change 
in  group  velocity  is  a  large  effect;  thus  it  should  be  readily  observable  in  experiments. 
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Figure  6.2:  a)  Two  photon  pulse  after  traveling  through  the  medium.  We  took  a 
gaussian  input  pulse  of  width  5(fl2/r)_1,  otherwise  the  parameters  are  as  in  Fig.  6.1b. 
b)  Comparison  between  numerical  simulations  of  the  change  in  group  velocity  vs  the 
change  in  group  velocity  for  the  NLSE.  They  agree  well  at  small  ODi „  but  quickly 
diverge  as  ODi  increases. 


6.3.2  Three  Photon  Solitons 

In  addition  to  the  two  photon  bound  states  there  should  also  be  a  manifold  of 
three  photon  bound  states  (Ben-Aryeh,  1999).  We  looked  for  these  by  doing  numerical 
simulations  of  the  steady  state  solution  for  three  photons.  The  results  are  shown  in 
Fig.  6.3  in  the  three  time  correlation  function  of  the  three  photons 


g(3\t1,t2,t3) 


(W(h)jv2(t2)iV3fe)) 

(iVi  (ti ) )  (iV2  (t2) )  (iV3  (t3) ) 


(6.15) 


where  IVj(tj)  are  photon  number  operators.  Fig  6.3a  shows  the  case  where  A  =  0 
where  we  expect  no  bound  state,  but  instead  the  three  photon  version  of  the  dissipa¬ 
tive  blockade  reported  by  Peyroncl  et  al.  (2012).  Fig.  6.3b  is  for  the  same  parameters 
as  Fig.  6.1b  with  A  =  2r  where  we  see  a  distinct  peak  when  all  three  photons  arrive 
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Figure  6.3:  a)  Three  time  correlation  function  g^\t2  —  ti,  t3  —  t\)  in  steady  state  with 
A  =  0.  The  lines  of  reduced  probability  correspond  to  regions  where  two  or  more 
photons  overlap  in  the  medium.  Parameters  are  as  in  Fig.  6.1b  except  with  A  =  0  b) 
Same  parameters  as  (a)  except  A  =  2T.  The  peak  in  the  middle  corresponds  to  the 
three  photon  bound  state  while  the  additional  features  arise  from  two  photon  bound 
states. 


simultaneously.  We  identify  this  peak  with  the  three  photon  bound  state. 


6.4  Conclusions 

We  have  studied  the  dynamics  and  formation  of  two  and  three  photon  solitons  in 
a  strongly  interacting  Rydberg-polariton  gas.  Such  states  may  be  useful  for  achieving 
quantum  gates  between  photons  and  creating  large  entangled  states.  In  future  work 
we  will  add  repulsion  interactions,  which  may  drive  the  system  to  a  crystalline  state 
of  photons. 
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Appendix  A 

Appendices  to  Chapter  2 

A.l  Parameters  Used  in  Simulations 

In  table  A.l  below  we  provide  a  summary  of  the  parameters  used  in  the  simulations 
for  each  figure.  While  many  parameters  are  chosen  to  be  consistent  with  experiments, 
not  all  those  presented  are  self-consistent  or  experimentally  realistic.  In  particular, 
in  Fig.  2.5e  the  A0  parameter  is  unphysically  large  and  in  Figures  2.4,  2.5ab  and  2.8 
the  small  mo  values  correspond  to  very  large  magnetic  fields. 

A. 2  Variables 

In  this  appendix  we  describe  a  systematic  approach  to  coarse  graining  the  electron 
wavefunction  in  solving  the  semiclassical  equations  of  motion,  which  we  refer  to  as 
the  Independent  Random  Variable  Annular  Approximation  (IRVAA).  We  construct 
a  sequence  of  discretizations  of  the  wavefunction  for  which  we  can  provide  a  rigorous 
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Table  A.l:  Parameters  used  in  the  simulations  shown  in  the  figures  of  Chapter  2  and 
Appendix  A. 


Fig. 

Ao 

r0 

A_ 

A+ 

Ao 

r* 

m0 

V 

M 

2.3 

0.5 

0.005  -  0.5 

O 

1 

o 

0 

0 

0 

0 

0.005 

400 

2.4 

5  •  mo 

/c/2 

1.25  •  mo 

5 

5  •  mo 

2.7 -To 

10"3  -  10"1 

0 

400 

2.5ab 

0.19 

1 

0.0048 

5.8 

0.002 

2.7 

5  •  10~4 

0 

100 

2.5c 

0.78 

1 

0.19 

5.8 

0.08 

1 

0.01 

0 

100 

2.5e 

1 

1 

0.25 

0.5 

1 

1 

0.05 

0 

100 

2.6a 

1.99 

1 

0.143 

626 

0.5 

2.7 

0.01 

0 

100 

2.6b 

1.99 

1 

0.143 

626 

0.5 

2.7 

0.01 

10-4 

100 

2.7a 

0.014 

0.36 

0.0034 

5 

0.014 

2.7 -To 

0.0027 

4  •  10"4 

200 

2.7b 

0.013 

2.1 

0.0034 

5 

0.013 

2.7 -To 

0.0027 

2  •  10~3 

200 

2.8 

5  •  mo 

/c/2 

1.25  •  mo 

5 

5  •  mo 

2.7 -To 

o 

1 

CO 

1 

O 

1 

10-4 

400 

A. la 

0 

1 

0 

0 

0 

0 

0 

10“3 

600 

A.  lb 

0.5 

0.005  -  0.5 

0-0.4 

0 

0 

0 

0 

5  •  10“5 

1200 

bound  on  the  error  in  time  evolution  compared  to  the  exact  solution.  In  the  process 
we  also  introduce  a  new  set  of  statistically  independent  nuclear  spin  variables,  which 
are  a  convenient  basis  for  numerical  simulations. 

We  see  from  Eqs.  2.4  and  2.6  that  the  semiclassical  evolution  of  each  spin  depends 
only  on  the  vectors  L  and  R  (or  equivalently  on  D  and  S).  That  is,  if  we  know  Pd(t) 
(which  depends  only  on  L  and  R),  then  we  can  solve  for  the  dynamics  of  the  entire 
system.  However,  even  if  we  know  Pd(t ),  if  we  look  at  the  equation  of  motion  for  L 
we  find  that  it  generates  an  infinite  hierarchy  of  equations 

§=P,xL-,  (A.1) 

where  we  defined  L*  =  J2kghhi-  Now  L*  couples  to  the  variable  J2k9kAki  and  so  on. 


97 


116 


Appendix  A:  Appendices  to  Chapter  2 

To  find  an  approximate  solution  to  the  dynamics  we  would  like  to  find  an  effective 
method  to  truncate  this  infinite  hierarchy  of  equations.  For  simplicity  we  focus  on 
the  case  where  Pi  is  only  a  function  of  L,  reducing  it  to  a  single  dot  problem,  and 
drop  the  dot  indices  in  the  following  discussion.  We  also  work  in  the  continuum  limit, 
which  is  defined  by  a  nuclear  angular  momentum  density  I(r,  t )  =  i k(t)S(r  —  r^). 

Each  variable  in  the  hierarchy  of  equations  of  motion  (as  in  Eq.  A.l)  can  be 
expressed  as  an  integral 


®{t)  =  I  ddrg{r)ip(g{r))I{r,t\  (A.2) 

where  <p(x)  is  a  polynomial  in  x.  That  is,  there  is  a  one-to-one  correspondence 
between  polynomials  <p(x)  and  the  variables  in  the  EOM.  For  example,  L  corresponds 
to  4>{x)  =  1. 

We  would  like  to  think  of  a  truncation  procedure  as  any  procedure  that  provides 
a  reduced,  self-consistent  set  of  equations  describing  the  evolution  of  P,  equivalently 
L.  We  make  a  formal  definition  of  a  truncation  procedure  as  a  procedure  producing 
a  set  of  variables  k  —  1, . . . ,  M,  of  the  form  above  and  an  M  x  M  matrix  Q,  such 
that  =  L  and 

=  2^  Px  Qki^e- 

t 

Since  we  always  constrain  =  L,  we  always  have  =  1. 

To  construct  a  convenient  basis  of  nuclear  spin  variables  we  first  define  a  norm  (•) 


based  on  the  statistical  average  of  a  nuclear  spin  variable  in  the  infinite  temperature 
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ensemble,  i.e. 


<4  ■*>«•=/ <f,r  ddr'  g2(r)f(g(r))ip(g(r'))  { J(r )  •  J(r'))e 

=  J  d*r g2(rMg(r))i,(g(r))  (A.3) 

where  a  is  the  lattice  spacing,  (-)e  is  the  ensemble  average  over  the  initial  thermal  state 
and  we  took  (I(r)  ■  I(r'))  =  1(1 +l)8(r  —  r')  /  ad .  Now  we  can  construct  an  orthogonal 
set  of  polynomials  with  respect  to  this  norm  by  using  the  standard  Gram-Schmidt  pro¬ 
cedure  starting  from  the  polynomial  1.  This  gives  a  set  of  orthogonal  polynomials  ip>k 
and  associated  nuclear  spin  variables  <f>k  =  f  ddr  g(r)(pk(g(r))I(r,  t ),  which  are  statis¬ 
tically  independent  in  the  infinite  temperature  ensemble  (i.e.,  (<&/,.  •  <f>/)  =  3f^<5fcz)and 
satisfy  <f>i  =  L. 

The  equations  of  motion  (EOM)  for  these  variables  can  be  written  as 

P  X  Qnm^m  (A. 4) 

where  the  matrix  Qmn  is  a  tridiagonal  matrix  defined  by  the  recurrence  relations 

x<pn(x)  =  Q 

nn—l^Pn—1  Qnn^Pn  +  Qnn+l^n+l  (A. 5) 

and  we  used  the  fact  that  xipn(x)  only  has  a  non-zero  overlap  with  ipn  and  (pn±i. 

We  now  define  an  Mth  order  truncation  procedure  with  respect  to  the  variables 
<Eq  by  setting  Qmm+i  =  0.  The  central  result  of  this  appendix  is  encapsulated  by  the 
following  theorem  for  this  truncation  procedure. 
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Theorem:  For  a  given  wavefunction  g(r)  and  £  >  0,  the  above  truncation 

procedure  at  order  M  will  produce  an  effective  LM(t)such  that  \ L(t)  —  LM(t)\  <  e 
for  all  t  <  tM,  where  tM  is  a  time  scale  that  increases  linearly  with  M  and  L(t)  is 
the  exact  result  for  the  untruncated  system. 

We  begin  our  analysis  by  proving  that  any  truncation  procedure  is  equivalent  to 
a  discretization  of  the  function  g(r )  (i.e.,  an  annular  approximation),  by  which  we 
mean  a  representation  of  L  as 

M 

L  =  '£/g(rk)ik,  (A. 6) 

fc=i 

where  Ik  is  a  rescaled  nuclear  spin  variable  associated  with  position  rk. 

The  reverse  implication  is  clear  because  if  we  start  with  such  a  discrete  represen¬ 
tation,  then  the  variable  associated  with  the  polynomial 

M 

w(x )  =  ~g(rk)\ 

k= 1 

is  identically  zero.  That  is,  if  there  are  only  M  discrete  spins  in  the  system,  then 
there  are  only  M  statistically  independent  variables  &k  in  the  system,  and  <E>m+i  is 
naturally  zero.  This  result  naturally  truncates  Eq.  A. 4.  Consequently,  if  we  consider 
any  basis  of  polynomials  of  degree  less  than  M  and  its  associated  set  of  spin  variables, 
then  we  can  obtain  a  finite,  self-consistent  set  of  equations  for  the  evolution  of  L. 

The  forward  implication  follows  along  similar  lines.  If  M  —  1  is  the  maximal  degree 
of  the  set  of  polynomials  {(pk(x)}  associated  with  the  truncation  variables  {*£&}  and 
<f>M  is  the  spin  variable  corresponding  to  this  polynomial,  then,  when  we  compare 
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to  the  continuum  limit,  we  find  that  the  statement  that  d$M/dt  does  not  couple  to 
higher  degree  polynomial  variables  implies  the  existence  of  a  degree-M  polynomial 
w(x )  such  that 

J  ddr  g{r)  w(g(r))  I(r,  t)  =  0, 

for  any  J(r,f).  The  existence  of  such  a  polynomial  immediately  implies  that  we  can 
represent  L  in  the  discretized  form  of  Eq.  A. 6. 

We  have  now  reduced  the  problem  of  finding  an  optimal  truncation  procedure  to 
the  problem  of  finding  an  optimal  discretization  procedure  for  integrals  of  the  form 


j  ddr  g(r)ip(g(r))  I(r,t), 

where  ip(x )  is  a  polynomial  in  x.  Fortunately,  this  last  problem  is  solved  through 
the  theory  of  Gaussian  quadrature.  (Kress,  1998)  First,  though,  we  assume  that  our 
function  g(r)  is  spherically  symmetric  so  that  we  can  write  our  integrals  as  effective 
one-dimensional  integrals  with  respect  to  the  rescaled  angular  momentum  density 


J(r,  t) 


j  dQ  ad~1N(r)I{r,n,t)/S(d ) 


(A.7) 


where  f2  parameterizes  the  surface  of  a  d-dimensional  sphere,  a  is  the  lattice  spacing, 
S(d)  is  the  surface  area  of  a  unit  sphere  in  d  dimensions,  and  N{r )  =  S(d)  rd~1/ad-1  is 
the  number  of  nuclear  spins  at  radius  r;  for  example,  in  two  dimensions  N(r)  =  27ir/a. 
The  ensemble  average  of  J(r,f)  is  given  by  (J(r)  •  I{r'))  =  /(/  +  1  )N{r)8{r  —  r')/a. 
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To  begin  constructing  our  Gaussian  quadrature  rules  we  rewrite 


drN(r)g2{r)(p{g(r )) 


I(r,  t) 

N(r)g(r ) 


dx  u(x)(p(x) 


1(9  1(x),t) 

N(^g~1(x))x 


(A.8) 


where  x  =  g{r)  and  uj(x)  =  ^|9-i (g~1(x))x2  is  the  weight  function.  Standard 
results  in  the  theory  of  numerical  integration  imply  the  existence  of  a  set  of  orthogonal 
polynomials,  ipn,  with  respect  to  the  inner  product 


(M) 


dxuj(x )  f(x)  h(x) 


(A.9) 


such  that,  for  any  function  f(x),  the  Mth  order  quadrature  approximation  is  given 
by 

,1  m 

/  dxu(x)  f(x)  &  y ^ukf(xk),  (A. 10) 

Jo  k=  1 

where  Xk  are  the  zeros  of  p>M  and  the  weights  04  are  determined  by  the  condition  that 
Eq.  A.  10  is  exact  for  all  polynomials  of  degree  strictly  less  than  2 M.  The  error  in  this 
formula  decreases  exponentially  in  M,  or  better,  provided  that  /  is  smooth.  (Kress, 
1998)  In  addition,  these  polynomials  are  exactly  the  ones  we  used  to  construct  our 
truncation  procedure.  Consequently,  our  truncation  procedure  defined  above  is  equiv¬ 
alent  to  approximating  L  in  quadrature  as  in  Eq.  A. 6  with  Ik  =  Wkl(rk,t) / gkN (rk) ■ 
To  prove  the  theorem  we  first  note  that  from  the  definition  \P(L)\  <  1  for  all 
L.  Now  let  p  >  0  be  such  that  | P(L)  —  P(L')\<  p\L  —  L  \  for  all  L  and  L  .  We 
define  Ln(t)  =  f  ddrgn(r)I(r,t )  and  L^(t)  is  the  solution  for  the  equivalent  variable 
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in  the  truncated  system  of  equations.  To  provide  bounds  on  the  error  propagation 
we  define  8^f(t)  =  \Ln(t )  —  L^(t)\.  We  work  in  time  units  where  ma xrg[r)  =  1  and 
let  b  =  maxnj\Ln(t)\<  f  ddrg(r)(I  +  1).  Now  it  is  straightforward  to  show  that 

%  <  pbS r  +  (1  +  AM)Ci  <  C(5f  +  Cl)  (A.ll) 

where  (  =  max(pb,  1  +  pe)  and,  by  assumption,  we  are  restricted  to  short  enough 
times  that  8™  <  £.  By  construction,  8^(0)  —  0  for  n  <  M  while  for  n  >  M  8 „  is 
bounded  by  the  quadrature  error  on  the  integral  f  ddrgn(r)I(r,0),  which  is  less  than 
ce~M  for  a  constant  c  independent  of  M.  Using  Eq.  A.ll  we  can  then  bound  the 
error  on  8 <  ce_M(e2^  —  1).  This  implies  that  the  time  to  make  an  error  of  size 
£  scales  as  (1/2^)  log (eeM/c+  1)  rs_/  (M  —  logc/e)/2(t  for  large  M.  This  proves  the 
theorem. 

For  the  two  dimensional  Gaussian  g(r)  oc  e~r~^2a~  the  weight  function  w(x )  =  x 
and  the  associated  orthogonal  polynomials  are  the  Jacobi  polynomials.  The  matrix 
Q  is  then  given  by  standard  recurrence  relations  for  Jacobi  polynomials.  Once  the 
recurrence  relations  are  known,  one  can  work  with  the  <h- variables  without  converting 
between  the  original  nuclear  spin  variables  because  the  $  variables  were  defined  such 
that  they  are  initially  statistically  independent.  This  is  a  convenient  numerical  ap¬ 
proach  for  these  types  of  central  spin  problems,  and  it  was  used  in  all  of  the  numerics 
in  this  work. 
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Table  A. 2:  Relative  population  of  the  nuclear  species  xa,  effective  hyperfine  field  clue 
to  species  a  ba,  and  the  gyromagnetic  ratio  ya,  for  the  three  nuclear  species  in  GaAs. 


75  As 

69  Ga 

71  Ga 

xa 

1 

0.6 

0.4 

ba  (T) 

-1.84 

-1.52 

-1.95 

ry 

'a  ^  mT  ) 

45.96 

64.39 

81.81 

A. 3  Multiple  Nuclear  Species 

In  this  appendix  we  include  the  effects  of  multiple  nuclear  species  in  our  simula¬ 
tions  and  find  that  the  main  results  for  both  asymmetric  and  identical  results  carry 
through  much  the  same.  First  we  show  how  to  include  multiple  species  in  terms  of 
the  collective  <h-variables  and  then  we  present  the  simulation  results. 

When  multiple  species  are  taken  into  account  we  must  include  the  Larmor  pre¬ 
cession  of  the  nuclear  spins.  In  this  case  the  EOM  take  the  form 

ikd  =  Aeba  V0  l^kdl1  Pd  xl%d-uazx  I£d,  (A.  12) 

where  a  is  a  species  index,  u>a  =  7 a  BextT / ra  is  the  effective  Larmor  frequency,  ba  is 
the  bare  hyperfine  field  of  species  a,  7 Q  is  the  gyromagnetic  ratio  of  species  a,  Bext 
is  the  external  magnetic  field,  and  we  have  explicitly  included  the  factor  T / t0  ,  where 
T  is  the  total  time  of  the  nuclear  pump  cycle  and  ra  is  the  adiabatic  sweep  time. 

We  introduce  the  projector  function  n%d,  such  that  n%d  =  1  if  there  is  species  a  in 
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|  Q 

T 
qn 


0  0.2  0.4  0.6  0.S 

A_/A0 


Figure  A.l:  a)  As  in  Fig.  2.7,  with  parameters  chosen  as  in  Fig.  2  of  Gullans  et  al. 
(2010),  except  with  three  species.  Due  to  the  computational  cost  of  running  three 
species  of  spins,  simulations  were  run  for  only  10%  as  long,  and  the  range  of  Dz/Sz  is 
larger  as  a  result.  The  trend  that  Dz/Sz  is  in  good  agreement  with  the  single-species 
prediction  is  clearly  visible,  b)  Phase  diagram  with  multiple  species  and  mo  =  0. 


unit  cell  k  and  0  otherwise.  This  allows  us  to  write 


L  =  Y1  7ebav0  KM2  K(  Jl 

k,ct 


Q 

k£ 


b«  9u  Ike 


(A. 13) 


Here  we  have  defined  f b  to  be  the  standard  deviation  of  in  the  infinite  temperatme 

state,  exphcitly 


{Ikd  ■  Ik'd')  =  I  {I  +  l)4fc'<W<W,  (A-14) 

°2  =  < L 2)  /3  =  ^  7 ?b2a  xa  vl  iV'fcfl4  ~  (A.15) 

A:, a 

where  xQ  =  (7r£d)  ts  the  relative  proportion  of  species  a  on  the  sites  it  can  occupy, 
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9kd  oc  Vq  \^kd\2  are  chosen  to  satisfy  J2kglI(I  +  1)  =  3,  and  /  is  the  total  spin  of  a 
single  nuclear  spin  (/  =  3/2  for  all  species  in  GaAs). 

We  define  the  variables 


$ 


a 

n 


7T Ike, 

k 


(A. 16) 


where  tpen(x)  are  defined  as  in  Appendix  A. 2  and  are  independent  of  the  species,  i.e. 
lPq(x)  =  1  and 


^gldViigkd)  <pem(gkd)I{I  +  1)  =  3  8nm.  (A. 17) 

k 

These  definitions  have  the  implication  that  •  L")  ,)  =  <5nnA/t/t Aaa' ,  and  we  can 
draw  initial  values  for  each  of  them  from  a  normal  distribution.  Furthermore,  we  can 
express 


L  = 


Vie 


Y,  hay/x^ 


VT,abl 

All  these  definitions  are  equivalent  for  the  right  dot. 
In  these  variables  the  EOM  take  the  form 


(A.18) 


<E»"  = 


7  eba 
N 


P(  X  Un^n-l  +  an®r 


+  £n+l<&n+l)  ~UaZX 


(A. 19) 


where  we  have  used  the  definition  N_1  =  max*.  vq  \^kd\2  to  represent  the  number  of 
nuclear  spins  with  which  the  electron  has  significant  overlap.  For  a  two  dimensional 
gaussian  wave  function  we  have  N  =  2/3  'YhaxaAeb2aI(I  +  l)/fl2 

In  Fig.  A.l  we  include  the  three  nuclear  species  in  the  simulation  and  show  that 
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qualitatively  the  results  from  the  single  species  case  still  hold.  Fig.  A. la  shows  the 
asymptotic  ratio  of  Dz/Sz  as  the  relative  dot  sizes  are  varied,  where  we  see  good 
agreement  with  the  simple  prediction  given  in  the  introduction.  In  Fig.  A.  lb  we 
extract  the  phase  diagram  in  the  simplified  model  with  only  A0,_  and  r0  non-zero,  as 
in  the  model  of  Gullans  et  al.  (2010).  As  in  the  single-spin  case,  we  find  a  saturation 
regime  at  high  values  of  To/Ao  and  an  instability  regime  at  lower  values.  Unlike  in  the 
single-spin  case,  the  saturation  regime  does  not  broaden  at  higher  values  of  A_/A0. 
The  dashed  line  is  the  same  as  that  in  Fig.  2.3,  showing  the  simple  prediction  for  the 
phase  boundary  with  a  single  species,  from  (Gullans  et  al.,  2010).  The  lower-left  side 
of  the  phase  diagram  (the  region  most  easily  reached  in  experiments)  is  well-described 
by  this  prediction,  even  with  multiple  species. 
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B.l  Van  der  Waals  Interaction  with  the  Nanosphere 

A  ground  state  atom  experiences  an  attractive  van  der  Waals  (vdw)  force  when 
placed  near  the  sphere  due  to  the  virtual  emission  and  reabsorption  of  photons  re¬ 
flected  from  the  surface  (Wylie  and  Sipe,  1984).  This  is  a  purely  quantum  mechanical 
effect  and  can  be  interpreted  as  a  modification  of  the  Lamb  shift  due  to  the  presence 
of  the  material,  which  changes  the  photon  density  of  states.  In  particular,  if  we  write 
the  atom-photon  interaction  Hamiltonian  as 

Hj  =  -ix  •  E{r0)  (B.l) 


108 


127 


Appendix  B:  Appendics  to  Chapter  3 


where  /i  is  the  dipole  operator  and  E  is  the  electric  held,  then  using  second  order 
perturbation  theory  one  can  write  the  energy  shift  of  the  ground  state  as 


6Ea 


1  y  (0 j E -a \k)  ik I  1°)  (ffl  ha  k)  (el  HP  | g) 

^  “  UJk+  UJe 

k.e 


(B-2) 


where  |0)  refers  to  the  vacuum,  | k)  to  a  one-photon  state  in  the  kt\i  mode  of  the 
system,  and  | g,  e)  are  the  ground  and  excited  states  of  the  atom.  Applying  Kramers- 
Kronig  relations  one  can  rewrite  this  as  (Wylie  and  Sipe,  1984) 


h  h 

5Ea  =  ~— Im  /  dojGafi(ro,r0-,u)aap(u)  =  /  dh,Ga/3(r0lr0]i^)aap(i^) 

471  J 0  27T  J o 

(B.3) 

We  have  defined  the  correlation  functions  for  the  electric  held  and  atomic  dipole 
moments 


Gafi(r,r'-,t)  =  i([Eol(r,t),Efl(r',0)])G(t)/h  (B.4) 

aap(t)  =  i  (g\  0)]  I g)  Q(t)/h  (B.5) 

with  0(f)  is  the  Heaviside  step  function.  These  can  be  identified  with  the  held  and 
atomic  susceptibilities,  respectively.  The  held  susceptibility  can  be  obtained  from 
the  classical  solution  for  the  electric  held  of  an  oscillating  dipole  near  the  sphere 
(Wylie  and  Sipe,  1984).  The  van  der  Waals  interaction  is  obtained  from  the  rehected 
contribution  to  Gap.  We  work  in  the  quasistatic  limit  where  the  distance  between  the 
atom  and  sphere  is  much  less  than  a  wavelength.  This  results  in  the  rehected  held  of 
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a  dipole  p  above  sphere:  Er(r,r'-,u>)  =  —  V(p  ■  V')<hr(r,  r'\ u>),  where 


$r(r,r';w)  = 


47T6 


E 


e(co)  -  1 


a 


2n+l 


o  —  s(u)  +  1  +  1  jn  r,n+1  rn+1 


P, 


cos 


(«  -  »')) 


(B.6) 


£  is  the  dielectric  constant  of  the  sphere,  a  is  the  radius,  r’  is  the  position  of  the 
dipole,  and  Pn  is  the  nth  order  Legendre  polynomial.  The  reflected  greens  function 
is  defined  by  the  relation 


E?(r,r'-,u)  =  Graf}(r,r';uj)-pf)  (B.7) 

Gap{r,  r'\ u)  =  -VaV'/3®r(r,  r';  uj)  (B.8) 


The  van  der  Waals  force  for  a  ground  state  atom  is  dominated  by  the  exchange  of 
low-frequency,  off-resonant  photons.  This  is  to  be  contrasted  from  situation  for  the 
excited  states,  where  the  atom  can  emit  and  reabsorb  real  photons  at  the  resonance 
frequency  leading  to  an  additional  correction  to  the  van  der  Waals  force  (Chance 
et  ah,  1975).  Because  of  this  we  are  justified  in  taking  £  — y  —  oo  in  Eq.  B.6,  which 
allows  us  to  write 


Grzz(r,r) 

Grxx(r,r) 


Uv  aw 


Os 


1  a3  4  —  3 a2 /r2  +  a4/r4 
47re0  r6  (1  —  a2/r2)3 


yyV  ’  y  47re0r6(l-a2/r2)3’ 
C3  2a3 (6  —  3 (a/r)2  +  (a/r)4) 
r6  (l  —  (a/r)2)3 


h 


167T2e 


o  Jo 


d£a(i£) 


12  ' 


(B.9) 

(B.10) 

ft  7  2a6(6-3(a/r)2  +  (a/r)4) 

16  ft3a3  r6(l  -  (a/r)2)3 

(B.ll) 

(B.12) 
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where  7  is  the  spontaneous  emission  rate  for  the  two-level  atom  in  free  space.  In  the 
limit  a  <C  r,  C/Vdw  ~  1/r6,  as  expected  because  the  sphere  responds  like  a  dipole.  In 
the  opposite  limit,  when  (r  -  o)  «  a  we  reproduce  the  well  known  formula  for  the 
ground  state  shift  of  an  atom  above  a  perfectly  conducting  plane  t/vdw  =  C3/(r  —  a)3. 
For  Rb87,  A  780  nm  and  7  =  6  MHz,  if  we  take  a  sphere  with  a  20  nm  radius  this 
gives  the  typical  scale  for  Uv<m  ~  100  MHz,  which  is  quite  substantial. 

B.2  Heating  Rate  from  Inelastic  Light  Scattering 

Here  we  calculate  the  heating  rate  due  the  inelastic  light  scattering  from  the 
trapping  laser  including  the  interaction  with  the  nanosphere.  Because  of  the  tight  trap 
confinement  the  change  in  motional  state  arises  from  events  where  a  single  phonon  is 
added  or  subtracted  to  the  system  (Grimm  et  ah,  2000).  Expanding  the  fields  around 
the  trap  center  gives  the  heating  rate 

1  jump  -L  tot  co  ?  113.10  J 

TVjJT,z  0Z 

where  P  is  the  Rabi  frequency  of  the  trapping  light,  5  is  the  trapping  laser  detuning 
from  the  atomic  resonance,  c ot,z  is  the  trap  frequency,  E'R  =  9h2/2mzR  is  an  enhanced 
recoil  energy  due  to  the  tight  trap,  m  is  the  mass  of  the  atom,  and  rtot  is  the  total 
spontaneous  emission  rate  of  the  atom  including  both  radiative  emission  and  non- 
radiative  emission  into  the  surface  plasmon  modes  of  the  sphere.  The  lifetime  of  the 
trap  is  approximately  given  by  the  time  it  takes  for  the  atom  to  hop  out  of  the  trap 
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clue  to  such  absorption  processes 


nn2/6 

Tjump  hi^T,Z 


(B.14) 


We  express  rtot  =  rrad  +  rnon_rad  in  terms  of  both  radiative  and  non-radiative 
contributions.  The  radiative  contribution  can  be  found  from  the  dipole  moment 
induced  in  the  sphere  from  the  excited  atom 


r  rad  —  'y 


ft  + 


47re0^3 


(3(/t  •  z)  z 


2 


A) 


(B.15) 


The  non-radiative  emission  arises  from  near  field  coupling  of  the  atom  to  plasmon 
modes  of  the  sphere.  It  can  be  expressed  as  Tnon_rad  oc  Im(jp  ■  Er(r',r')),  where  Er 
is  the  field  calculated  in  Eq.  B.6.  Tnon_rad  contains  both  a  resonant  and  non-resonant 
contributions  from  the  dipole  and  multipole  contributions,  respectively 

Tnon-rad  _  6  a3  Im(a(u;a))  3  1  /e  -  1\  a8  (9  -  ll(a/r)2  +  4(a/r)4) 

T0  “fc3r6  4vr e0a3  +  2  £;3a3  m \e  +  1 )  C  (l  -  (a/r)2)3  J 

For  moderate  distances  from  the  sphere  we  see  that  Tnon.rad  is  dominated  by  the 
emission  into  the  resonant  surface  plasmon  mode.  In  addition,  this  emission  can  be 
substantially  greater  than  the  radiative  emission. 


B.3  Tuning  the  Lattice  Potential 

In  order  to  control  the  tunneling  rate  in  the  Hubbard  model,  one  needs  control 
over  the  trapping  potential  in  the  plane  of  the  lattice.  This  can  be  achieved  through 
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Figure  B.l:  a-b)  Contours  of  atomic  potential  in  MHz  for  a  ID  chain  of  silver 
nanoshells  including  vdw  with  light  blue-detuned  to  the  plasmon  resonance,  linear 
polarized  light  is  applied  from  the  side  and  circularly  polarized  light  is  applied  from 
above.  The  lattice  potential  can  be  tuned  by  changing  the  polarization  between 
linear  and  circular:  UZ/UQ  =  1  in  (a),  while  Uz/U0  =  0.75  in  (b).  c)  Lattice  poten¬ 
tial  along  the  chain  for  different  amounts  of  circular  polariztion:  Uz/U0  =  1  (solid), 
Uz/U0  =  0.75  (dashed),  and  Uz/U0  =  0.5  (grey). 


polarization  control  similarly  to  the  loading  procedure.  Figure  B.4  demonstrates  this 
tuning  in  a  lattice  formed  by  linearly  polarized  light.  Here  adding  circularly  polarized 
light  lowers  the  potential  in  the  plane  of  the  lattice,  while  simultaneously  maintaining 
the  trap  in  the  vertical  direction. 


B.4  Effective  Scattering  Length  in  Tight  Traps 

The  scattering  problem  for  two  atoms  in  a  three-dimensional  isotropic  trap  in¬ 
teracting  via  a  contact  potential  can  be  solved  exactly.  We  follow  the  approach  as 
described  by  Busch  et  al.  (1998)  and  Bolda  et  al.  (2002)  and  define  an  energy  depen¬ 
dent  effective  scattering  length  as  aeff(E)  =  -  tan  r]o(k) /k.  We  find  the  eigenvalues 
of  the  system  by  solving: 

— j —  =  HE)  (B.17) 
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where  l  =  \Jh/muj  is  the  harmonic  oscillator  groundstate  length  and  the  so  called 
’intercept’  function  f(E)  is  defined  as: 


(  IX E  7r\ 

\2hLj  +  AJ 


r(^  +  !) 


We  calculate  the  effective  scattering  length  by  using  the  accumulated  phase  method 
as  described  by  (Verhaar  et  ah,  2009);  we  solve  the  radial  Schrodinger  equation  be¬ 
tween  r  =  a%n  =  20  ao  and  r  — >•  oo  where  we  apply  the  known  scattering  length  as  a 
boundary  condition  at  r  — >•  oo,  this  gives  us  the  phase  of  the  wavefunction  at  r  =  ain. 
Subsequently  we  calculate  the  effective  scattering  length  as  a  function  of  energy  E  by 
using  the  phase  at  r  =  ajn  as  the  boundary  condition.  We  assume  the  accumulated 
phase  is  energy  independent  over  the  energy  range  we  consider.  This  results  in  an 
energy  dependent  scattering  length.  We  verified  the  validity  of  the  accumulated  phase 
method  by  comparing  to  the  results  for  23Na  obtained  by  Bolda,  et.  al.  and  find  good 
agreement  (Bolda  et  ah,  2002).  The  approach  breaks  down  if  the  harmonic  oscillator 
length  becomes  smaller  than  the  van  der  Waals  range  (/  <  rvcnv )  which  is  defined  as 
Tvdw  =  \  {^ixCe/h2)1^4 .  For  87Rb  this  implies  the  trapping  frequency  should  be  less 
than  12  MHz. 


Figure  B.4  shows  the  results  of  this  calculation  for  87Rb  with  a  1  MHz  trapping 
frequency.  We  took  a  triplet  scattering  length  of  =  98.99a0  and  C6  =  4698a0 
(van  Kempen  et  al.,  2002).  For  these  parameters  we  find  a  resonance  in  the  effective 
scattering  length  near  E  ~  h  x  9.5  MHz  ~  kBx450  pK,  which  is  between  the  4th  and 
the  5th  vibrational  state.  In  the  inset  to  Figure  3.2  we  show  the  effective  scattering 
length  for  the  lowest  vibrational  level  as  a  function  of  trap  frequency  where  we  see  a 
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Figure  B.2:  The  effective  scattering  length  (blue  solid  curve)  and  the  intercept  func¬ 
tion  (red  dashed  curve)  for  a  trap  frequency  of  u  =  1  MHz.  The  eigenvalues  of  this 
system  correspond  to  the  crossings  of  the  two  curves. 

resonance  at  u  ~  3.8  MHz. 

This  scattering  problem  will  be  also  affected  by  the  sphere  because  it  modifies  the 
vdw  interaction  between  the  atoms.  However,  the  spheres  contribution  will  be  small 
compared  the  bare  vdw,  provided  the  typical  distance  between  the  atoms  on  a  single 
site  is  much  less  than  their  distance  to  the  sphere. 

B.5  Two  Atom  Entanglement  on  the  Lattice 

For  two  atoms  on  sites  0  and  n  we  take  the  density  matrix  evolution 

H  =  +  O  +  -  alg  +  h  c- )  +  B^°ls  +  a< )s  +  h  c •)  (B-18) 

P  =  -*  \H-  P]  ~  70 nV[o]e  +  <Xge]p  ~  <57n  (Vlale}  +  V[(j]e\  +  V[(j\e]  +  V[(T*e\)p  (B.19) 
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Figure  B.3:  a)  Level  diagram  for  two  atoms  showing  transitions  driven  by  external 
fields  and  decay  pathways,  b)  Level  diagram  showing  effective  transition  rates  in  the 
ground  state  manifold.  The  pumping  rate  into  the  state  \sg)  —  |ps)  is  much  larger 
than  the  rate  out  of  it.  c)  Shows  the  infidelity  for  preparing  the  singlet  state  after 
optimizing  Bx  and  il  as  the  sub-radiant  states  linewidth  becomes  narrower. 


where  A  is  the  detuning  between  the  control  fields  and  the  excited  state,  Q  is  an 
optical  control  field,  Bx  is  a  transverse  ground  state  magnetic  field,  and  D[c\p  — 
1/2 {c'c.p)  —  cpC .  In  addition  to  the  decay  from  |e)  to  | g)  through  the  plasmons  we 
assume  there  is  an  additional  decay  from  \e)  to  |s)  that  occurs  at  the  rate  8jn.  This 
term  is  essential  to  remove  entropy  from  the  system  to  cool  into  the  singlet  state. 
The  relevant  process  are  shown  schematically  in  Fig.  B.5a. 

The  minimal  error  in  preparing  the  singlet  state  decreases  linearly  with  the  ratio 
djn/jnn  as  shown  in  Fig.  B.5c.  This  can  be  understood  in  the  limit  of  weak  driving 
Q  <C  8jn  7 nn-  In  this  limit  the  excited  states  can  be  adiabatically  eliminated  to 
give  the  effective  evolution  depicted  in  Fig.  B.5b.  Because  the  optical  pumping  rate 
out  of  a  state  increases  inversely  with  the  linewidth,  the  pumping  rate  into  the  singlet 
state  Ri„  «  CL2 /5^fn  can  be  much  larger  than  the  pumping  rate  out  of  it  Rout  ~  Q2/7„„. 
If,  in  addition,  Bx  2>  C2 /8^n  the  triplet  states  are  completely  mixed  and  all  triplet 
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states  can  be  optically  pumped  into  the  singlet  state.  The  ratio  Rout/Rm  rs_/  Sin/ Ann 
then  determines  the  relative  population  in  the  triplets  to  the  singlet  state,  giving  the 
hdclity  F  «  1  —  Sin/ Ann- 

As  a  remark  we  note  that  Eqs.  16  and  17  can  be  mapped  exactly  to  a  cavity  QED 
model  by  replacing  7 0n  by  g2  j  k  and  <5yn  by  7,  where  g  is  the  coupling  of  a  single 
atom  to  a  single  cavity  photon,  k  is  the  cavity  decay  rate,  and  7  is  the  free  space 
decay  rate.  In  this  case  the  fidelity  scales  as  1  —  1/P  where  P  =  g2  /  n 7  is  the  Purcell 
factor.  This  linear  scaling  of  the  singlet  fidelity  with  the  Purcell  factor  agrees  with 
the  limit  obtained  by  Kastoryano  et  al.  (2011)  using  a  similar  dissipative  approach. 
The  main  difference  between  the  two  schemes  is  that  for  the  scheme  by  Kastoryano 
et  al.  (2011),  the  cavity  resonance  is  assumed  to  be  far  detuned  from  the  atomic 
resonance,  while  in  the  present  approach  the  two  resonances  are  the  same.  Thus  they 
operate  in  qualitatively  different  regimes  of  cavity  QED.  In  the  off  resonance  case  the 
cavity  interaction  shifts  the  excited  state  energies  for  the  states  | eg)  ±  | ge),  while  in 
the  resonant  case  the  cavity  interaction  results  in  different  linewidths  for  | eg)  ±  | ge). 
Clearly  either  phenomenon  is  sufficient  for  ground  state  entanglement  generation. 
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C.l  Nonlinear  Conductivity 

The  nonlinearity  can  be  derived  from  the  Boltzmann  equation  for  the  2D  elec¬ 
tron  distribution  function  f(x,k,t )  and  Poisson’s  equation  for  the  electric  potential 


dtf  +  vFk-dxf  +  edxp  ■  dkf  =  0  (C.l) 

(dl  +  d%)<p  =  en5(z)/e0£  (C.2) 

where  n  =  j  dkf  is  the  2D  electron  density  dx  =  dxx  +  dyy  and  d k  =  d^xx  +  d kyy. 
Taking  x  to  be  the  propagation  direction  the  nonlinear  equations  for  the  moments  n 
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and  n  v  —  f  dkvpk  f  can  be  derived  as 


dtn  +  dxn  v  —  0 

dtv - -dx(p  +  vdxv+  — — — dx(p  Sn  =  0 

m*  2m*  n0 


(C.3) 

(C.4) 


where  the  effective  mass  for  the  plasmon  excitations  is  m*  =  hkp/vp,  Sn  =  n  —  n$  and 
n0  =  kp/n  is  the  equilibrium  electron  density.  Linearizing  these  equations  around  n0 
and  v  =  0  gives  the  plasmon  dispersion  from  Eq.  1  (Fetter,  1973).  The  nonlinearity 
is  described  by  the  last  two  terms  into  Eq.  C.4,  where  the  second  term  oc  dcpSn 
arises  from  the  linear  band  structure  in  graphene  and  is  absent  for  electrons  with  a 
parabolic  dispersion.  To  find  the  nonlinear  conductivity  we  can  expand  v  and  n  in 
spatial  Fourier  components  and  solve  the  resulting  coupled  equations  combined  with 
Poisson’s  equation  (Denardo  and  Putterman,  1988).  This  allows  us  to  express  <r3(a;) 
through  the  identity  env  =  cr(u)E  +  a3(co)E3. 

To  solve  for  the  nonlinear  shift  in  the  cavity  we  use  the  boundary  condition  that 
v  ■  n  —  0.  This  allows  us  to  represent  v  =  x  vp  sin px  and  n  =  nq  cos  Qx  where 
q  =  mksp  for  some  integer  m.  Inserting  this  solution  into  Eq.  C.3-C.4  leads  to  coupled 
nonlinear  equations  for  nq  and  vq 


E 


sin/jxl  vp  —  ,op 
n0p 


1 

2 


n„  ]  —  -  y ^IpVpVq  sin (77  -  q)x 


p,q 


PVpVq 


3  ^ p  npnq 


2  p  Uq 


sin(p  +  q)x] 


1 


(C.5) 


cos px(hp  +  nopvp)  =  -  npvq  (p  —  q )  cos(p  —  q)x  —  (p  +  q)  cos (p  +  q)x 


p,q 


119 


138 


Appendix  C:  Appendices  to  Chapter  4 

where  up  is  given  by  Eq.  1.  These  equations  can  be  solved  in  perturbation  theory  to 
End  the  nonlinear  frequency  shift  of  the  plasmon  resonance  given  in  Eq.  4. 


C.2  Quantizing  the  Plasmon  Mode 


To  quantize  the  plasmon  mode  we  use  the  Hamiltonian  (Gervasoni  and  Arista, 
2003) 


H 


-  I  dxeSncp  +  - 


*  ;r,2 


2  _ 
Am* 
4  n0 


dxn0m*v 


J2  5nl  +  5iil) 

<?  q 


(C.6) 


where  A  =  n 2 /k2sp  is  the  area  of  the  sheet  and  we  used  the  relation  vq  =  —Snq/qno 
from  the  continuity  equation.  This  Hamiltonian  can  be  quantized  in  the  usual  way 
by  defining  Snq  =  ^~(aq  +  cdq)  for  bosonic  operators  aq  such  that  aq  =  —iujqaq  and 
7 q  =  ‘Iq^/iOqpjpl'nA.  This  leads  directly  to  Eq.  4. 


C.3  Coupling  between  Nanoribbon  and  Cavity 

To  calculate  the  coupling  between  the  cavity  and  the  proximal  nanoribbion  we 
use  the  electric  potential  of  the  nanoribbon  plasmons  acting  on  the  graphene  cavity 


<Pr{x)  = 


Wr 


Y,Jdx 


,e  nrk  cos  kx' 
I  x  +  d 


x ' 


E 


e  n. 


(C.7) 


47re0£  ^  k2(x  +  d )5 


120 


139 


Appendix  C:  Appendices  to  Chapter  4 


where  we  assumed  d  S>  W,  Xsp.  Inserting  this  into  .  Eq.  C.4  gives  the  coupling  between 
each  plasmon  mode  k  in  the  nanoribbon  with  the  plasmon  mode  q  of  the  cavity  as 


^ kq 


1  krF  Wr 

1 

J  kcFL 

ujrk  +  c oc  q2k  d3 

(C.8) 


where  L  is  the  length  of  the  nanoribbon,  kpC  is  the  Fermi  wavevector  and  uJkl  is  the 
dispersion  of  the  ribbon(r)  and  cavity (c).  Applying  Fermi’s  Golden  rule  gives  the 
decay  rate  of  the  cavity  mode  into  the  nanoribbon  plasmons  given  in  Eq.  8. 
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D.l  Numerical  Code  for  Two-Photon  Time  De¬ 
pendent  Dynamics 


Here  we  show  the  numerical  MATLAB  code  we  wrote  for  time  dependent  dynamics 
of  2  photons. 

°/o%7o  Evolves  two  photon-atom  wavefunction  in  time  including  Rydberg 
7o7o7o  interactions 
7o7o7»  Variables: 

7o7«7o  PhotIn_x  -  Input  photon  wavefunction  is  a  function  of  space  x  and  integer  time 
7o7«7o  WFla  -  WF  in  region  with  1  photon  in  the  medium  and  one  behind 
7o7«7o  pWFlab  -  photon  WF  at  boundary  of  region  with  1  phot  out  and  1  phot  at 
7o7«7o  start  of  the  medium 

7o7«7o  WF2  -  WF  inside  medium  composed  of  3  level  Rydberg  atoms 

7o7«7o  RydV  -  Flag  which  says  whether  or  not  there  are  interactions 

7o7«7o  Rydlnd  -  Pairs  of  coord  separ.  by  Ryd  blockade  for  square  potential 

7o7«7o  u2RightB  -  Boundary  WF  for  region  with  1  phot  out  and  1  phot  in 

7o7«7o  u2upperB  -  Boundary  WF  for  region  with  1  phot  out  and  1  phot  at  start 

7o7«7o  of  medium 

7o7«7o  Uevolxl  -  Evolution  operator  e~{-i  HatomPhot  delt} 

7o7«7o  Uevolx2  -  Tensor  product  of  Uevolxl  with  itself 
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for  m=l:Nt 

%  Update  WF  in  region  with  1  phot  in  medium  and  1  behind 
WFla(l , : )=circshift (WFla(l , : ) , [0 ,nshift] ) ; 

WFla(l , 1 :nshift)=PhotIn_x(0 ,m) ; 

WFla=Uevolxl*WFla; 

°/  °/ o/o/ 0/0/ 0/0/0/ 

/o  /o  /o  /o  /  o  /o  /o  /o  /o 

%  Record  output  at  boundary  of  region  with  1  phot  out  and  1  phot  at 
%  start  of  the  medium 
pWFlaB=circshift(pWFlaB,-l) ; 
pWF 1 aB ( end) =WF 1 a (1 , end ) ; 

°/  7  °/  7  yrrrr/ 

/o  /  o  /o  /o  /o  /o  /o  /o  /o  /o 

%  Include  second  photon  to  get  two  photon  WF 
WFlaB=PhotIn_x(0 ,m) . *WFla; 

70/0/77770/77 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

*/.  Update  two  photon  WF  with  both  photons  in  medium 
WF2(1 )=circshift(WF2(l ,:,:), [0  nshift  nshif t] ) ; 

WF2(2, : , : )=circshift(WF2(2, : , : ) , [0  nshift  0]); 

WF2 (3 , : , : ) =circshif t (WF2 (3 [0  nshift  0] ) ; 

WF2d:3d,  : )  =WFlaB ; 

WF2 (4 , : , : ) =squeeze (WF2 (2 , 5 ; 

WF2 (7 , : , : ) =squeeze (WF2 (3 ; 

WF2(1,  :  d)=WF2(ld,  :); 

*/.  etprod  is  an  external  function  for  multidimensional  tensor  products 
*/.  Can  be  found  online 

WF2=etprod(,al2),Uevolx2, JabJ ,WF2, ’bl2;) ; 

'/.  Include  Rydberg  interaction 
if  RydV~=0 

WF2 (9 , Rydlnd) =0 ; 

end 

°/°/ °/°/ °/°/ °/“/ °/°/ 

/ 0  / 0  /O  /o  /0  /O  /0  /O  /O  /o 

*/.  Extract  boundary  conditions  for  exiting  the  medium 

u2rightB=WF2 ( [1  4  7],:, end); 

°/°/ °/°/ 0/0/°/ °/°/ 

/ 0  /o  /O  /O  /0  /O  /0  /O  /O 

7.  Update  WF  in  region  with  1  phot  out  and  1  phot  in  medium 
WF2b (1 , : , : )=circshif t (WF2b(l [0  1  1]); 

WF2b(2, : , : )=circshift (WF2b(2 , : , :) ,  [0  0  1]); 

WF2b(3, : , : )=circshift(WF2b(3, : , :) , [0  01]); 

u2upperB=zeros (1 , Nb) ; 

u2upperB(l :min(m,Nb) )=f liplr (pWFlaB(l :min(m,Nb) ) ) . *PhotIn_x(0 ,m) ; 
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WF2b(2 ,  1 ,  : )=PhotIn_x(0,m) . *WFla(2, end) ; 

WF2b(3, 1 , : )=PhotIn_x(0,m) . *WFla(3, end) ; 

WF2b (1,1, : ) =u2upperB ; 

WF2b ( : , : , 1 ) =u2r ightB ; 

WF2b=etprod( ’ al2 ’ ,Uevolxl, ’ab’ ,WF2b, ^12’) ; 

"/  0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 

/o  /  o  /o  /o  /  o  /o  /o  /o  /o  /o  /o  /o  /o 

l  Save  ouptut  by  recording  values  for  WF2b(l , end, : ) 

end 


D.2  Numerical  Code  for  Three-Photon  Steady  State 
Solution 


Here  we  show  the  numerical  MATLAB  code  we  wrote  for  steady  state  solutions 
of  up  to  three  photons. 

°/„°/o7o  Finds  wavefunction  (WF)  in  steady  state  for  three  photons  incident  on 
HI  the  medium 
III  Variables: 

HI  WF2-5  -  WF  in  regions  2-5,  4B,  and  5B  defined  below 

HI  eitR  -  Ratio  of  amplitude  between  S  and  E  state  for  dark  state 

III  polariton 

HI  alpha  -  Initial  amplitude  of  coherent  state  input 
III  Natoms  -  Number  of  atoms 

HI  Nmeas  -  Number  of  grid  points  taken  past  the  medium  where  photons  are 
HI  measured 

HI  Rydlnd  -  Pairs  of  coord  separ.  by  Ryd  blockade  for  square  potential 


III  2  photons  in  1  behind 
III  [EEE ; EES ; ESE ; ESS] ; 

WF2=zeros (4 , Natoms , Natoms) ; 

III  3  photons  in 

HI  [EEE ;  EES ;  ESE ;  ESS ;  SEE ;  SES ;  SSE ;  SSS]  ; 
WF3= (zeros (8 , Natoms , Natoms , Natoms) ) ; 


124 


143 


Appendix  D:  Appendices  to  Chapter  6 


HI  2  photons  in  1  out 

HI  [EEE ; EES ; ESE ; ESS ; SEE ; SES ; SSE ; SSS] ; 

WF4= (zeros (8 , Natoms , Natoms , Nmeas) ) ; 

HI  1  photon  in  2  out 
III  [EEE; SEE]; 

WF5= (zeros (2 , Natoms , Nmeas , Nmeas) ) ; 

III  1  photons  in  1  out  1  behind 
HI  [EEE ;  EES ;  ESE ;  ESS]  ; 

WF4B=zeros (4, Natoms, Nmeas) ; 


HI  2  photons  out  1  behind 
III  EEE 

WF5B=zeros (4 , Nmeas , Nmeas) ; 


HI  Parameters  governing  solution  in  region  1  with  1  photon  in  and  2  behind 

eitR=-gc/Omega; 

alpha=l ; 

1111111  Region  2  1111111111 


HI  Initial  conditions 
WF2(1 , 1 , : )=alpha~3; 

WF2(1, : , l)=alpha~3; 

WF2(2 , 1 , : )=eitR*alpha~3 ; 

WF2(3, : , l)=eitR*alpha~3; 

WF2(4, 1,1) =eitR/2 . * (WF2 (2,1,1) +WF2 (3,1,1)); 

WF2 (4 , Rydlnd) =0 ; 

HI  Update  WF  in  region  2 
for  L=l: Natoms 
if  L>1 

HI  EEE  has  3  spatial  derivatives  so  we  use  value  from  all  three 
HI  one  behind,  similarly  EES  and  ESE  have  2  spatial  derivatives 
HI  and  ESS  has  one . 

WF2(l,L,2:end)=WF2(l,L-l,l:end-l)+delt.*  . . . 

etprod( ’ al2 ’ ,HamEffx2(l , : ) , 5 ab; ,WF2( : ,L-1 , 1 : end-1) , ^12’) ; 
WF2(2 ,L, 2 : end) =WF2 (2 , L-l , 2 : end)+delt . *  . . . 

etprod( ’al2’ ,HamEffx2(2, : ) , ’ ab5 ,WF2( : ,L-l,2:end) , ;bl2’) ; 

WF2 (2 , L, 1) =WF2 (3 , 1 , L) ; 
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end 

for  k=l:Natoms 
if  k>l 

WF2(3,L,k)=WF2(3 ,L,k-l)+delt . *  . . . 

etprod( ’  al2’ , HamEffx2(3, :) , ’ab’ ,WF2(: ,L,k-l) , ;bl2’) ; 

end 

if  RydInd(L,k)==0 

WF2 (4, L , k) =eitR . * (WF2 (2 ,L ,k) +WF2 (3 ,L ,k) ) /2 ; 

end 

end 

end 

inini  Region  4b  nnnnn 

HI  Use  WF2  to  get  boundary  condition  for  region  4  with  2  in  1  out. 
U1  Initial  conditions 
WF4B ( : , : ,1)=WF2(: , : .end) ; 

WF4B (1 , 1 , : )=alpha~3; 

WF4B (2,1, : ) =WF4B (2,1,1)  ; 
for  L=l:Natoms 
if  L>1 

WF4B (3 , 1 ,L)=WF4B(3 , 1 , L-l)  +delt.*  ... 

HamEf fx2(2 , : ) *squeeze(WF4B( : , 1 ,L-1) ) ; 

end 

WF4B(4 , 1 , L)=eitR . * (WF4B (2,1, L)+WF4B (3 , 1 ,L) ) /2 ; 

end 

HI  Full  WF 
for  L=2:Nmeas 


WF4B ( 1, 2 : end, L)=WF4B( 1,1: end-1, L-l)  . . . 

+  delt . *HamEffx2 (1 , : ) ^squeeze (WF4B( : , 1 : end-1 , L-l) ) ; 
WF4B(3, : ,L)=WF4B(3, : ,L-1)  ... 

+  delt . *HamEffx2 (3, : ) ^squeeze (WF4B( : , : ,L— 1) ) ; 
for  k=l:Natoms 
if  k>l 

WF4B(2,k,L)=WF4B(2,k-l,L)  +delt.*  ... 
HamEffx2(2, : ) ^squeeze (WF4B ( : ,k-l,L)) ; 

end 

WF4B(4,k,L)=eitR.*(WF4B(2,k,L)+WF4B(3,k,L))/2; 

end 


end 
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1111111  Region  5B  1111111111 


HI  Find  BC  for  region  5  with  2  out  and  1  in  from  region  4B 
7o%7o  Initial  conditions 
WF5B ( : ,1, : ) =WF4B ( : ,end, :) ; 

WF5B (1 , : ,1)=WF5B(1,1, :)  ; 

WF5B (2 , : , 1) =WF5B (3 , 1 , :)  ; 

WF5B (3 , : , 1) =WF5B (2 , 1 ,  :)  ; 

WF5B (4 , : , : )=eitR*(WF5B(2, : , :)+WF5B(3, : , :))/2; 


for  L=2:Nmeas 

WF5B ( 1 , L , L : end) =WF5B ( 1 , L-l , L-l : end-1) ; 
WF5B ( 1 , L : end , L) =WF5B ( 1 , L , L : end) ; 

end 


7o7.7o7o7.7o7o  Region  3  1111111111 


7.7.7o  Initial  conditions 
WF3(1 : 4, 1 , : , :)=WF2; 

WF3(2 , : ,1, :)=WF2(2, :  ,  :); 

WF3(3, : , : ,1)=WF2(2, :  ,  :); 

WF3 ( 1 , : , 1 , : ) =WF2 ( 1 , : , :); 

WF3(1 , : , : , 1) =WF2 (1 , : ,  :); 

WF3(5, : ,1, : ) =squeeze (WF3 (3,1, :,:)); 

WF3(5, : , : , 1) =squeeze (WF3 (3,1, :,:)); 

WF3(6 , : , 1 , : ) =squeeze (WF3 (4,1, :,:)); 

WF3(7, : , : ,l)=squeeze(WF3(4,l, :,:)); 

7.7.7.  SSS  has  no  spatial  derivatives  so  steady  state  solution  can  be  derived 
WF3(8, 1 , 1 , l)=eitR/3 . *sum(squeeze (WF3( [4  6  7], 1,1,1))); 

WF3 (8 , Rydlndl23) =0 ; 

WF3 (4 , Rydlnd23) =0 ; 

WF3 (7 , Rydlndl2) =0 ; 

WF3 (6 , Rydlndl3) =0 ; 


7.7.7.  Initial  conditions 
for  k=l:Natoms 

for  m=l:Natoms 

if  m>l  &&  Rydlndl2(l,k,m)==0 

WF3(7, 1 ,k,m)=WF3(7, 1 ,k,m-l) . . . 

+delt . *HamEf fx3(7 , : ) *squeeze(WF3( : , 1 ,k,m-l) ) ; 

end 
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if  k>l 

if  Rydlndl3(l ,k,m)==0 

WF3(6, 1 ,k,m)=WF3(6, 1 ,k-l ,m) . . . 

+delt . *HamEffx3(6, : ) ^squeeze (WF3( : , 1 ,k-l ,m) ) ; 

end 
if  m>l 

WF3(5, 1 ,k,m)=WF3(5, 1 ,k-l ,m-l) . . . 

+delt . *HamEffx3(5 , : ) ^squeeze (WF3( : , 1 ,k-l ,m-l) ) ; 

end 

end 

if  Rydlndl23(l ,k,m)  ==  0 

WF3(8 , 1 ,k,m)=eitR/3 . *sum(squeeze(WF3( [4  6  7] , 1 ,k,m) ) ) ; 

end 

end 

end 


m  Initial  conditions 


WF3 (2 , 

,  : , l)=squeeze(WF3(5, 1 ,  : 

,:)) 

WF3 (3 , 

, 1 , : ) =squeeze (WF3 (5,1, : 

,:)) 

WF3 (4 , 

, : , l)=squeeze(WF3(6, 1 , : 

,:)) 

WF3 (4 , 

, 1 , : ) =squeeze (WF3 (4 , : , : 

,D) 

WF3 (6 , 

, : , l)=squeeze(WF3(7, 1 , : 

,:)) 

WF3 (7 , 

, 1 , : ) =squeeze (WF3 (7,1, : 

,:)) 

WF3 (8 , 

, : , l)=squeeze(WF3(8, 1 , : 

,:)) 

WF3 (8 , 

, 1 , : ) =squeeze (WF3 (8,1, : 

,:)) 

V/X  Find  full  WF  taking  into  account  when  2  phot  are  within  a  blockade 
Ul  Could  be  made  much  faster  with  more  efficient  matrix  representation 
111  or  by  writing  in  c  and  using  mex 

for  L=2:Natoms 

for  k=2:Natoms 

for  m=2:Natoms 

WF3(2,L,k,m)=WF3(2,L-l,k-l,m)  . . . 

+delt . *HamEf fx3(2 , : ) *squeeze(WF3( : , L-l ,k-l ,m) ) ; 
if  RydIndl3(L,k,m)==0 

WF3(6,L,k,m)=WF3(6,L,k-l,m) . . . 

+delt . *HamEf fx3(6 , : ) ^squeeze (WF3( : ,L,k-l ,m) ) ; 

end 

WF3 ( 1 , L , k , m) =WF3 ( 1 , L- 1 , k- 1 , m- 1 ) . . . 

+delt . *HamEf fx3(l , : ) *squeeze(WF3( : ,L-1 ,k-l ,m-l) ) ; 
WF3(5,L,k,m)=WF3(5,L,k-l,m-l) . . . 
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+delt . *HamEf fx3(5 , : ) *squeeze(WF3( : ,L,k-l ,m-l) ) ; 
if  RydInd23(L,k,m)==0 

WF3(4,L,k,m)=WF3(4,L-l,k,m) . . . 

+delt . *HamEf fx3(4, : ) ^squeeze (WF3( : ,L-1 ,k,m) ) ; 

end 

WF3(3,L,k,m)=WF3(3,L-l,k,m-l)  . . . 

+delt . *HamEf fx3(3 , : ) *squeeze(WF3( : ,  L-l ,k,m-l) ) ; 
if  RydIndl2(L,k,m)==0 

WF3(7,L,k,m)=WF3(7,L,k,m-l) . . . 

+delt . *HamEffx3(7, : ) ^squeeze (WF3( : ,L ,k,m-l) ) ; 


end 

if  RydIndl23(L,k,m)  ==  0 

WF3(8,L,k,m)=eitR/3. *sum(squeeze (WF3 ( [4  6  7] ,L,k,m) ) ) ; 

end 

end 


end 


end 


1111111  Region  4  1111111111 


HI  Initial  conditions 
WF4(: , : , : ,1)=WF3(: , : , : ,end) ; 

WF4(1,1,:,:)=WF4B(1 
WF4 ( 1 , : , 1 , : ) =WF4B ( 1 , : , 

WF4 (2 , 1 , : , : ) =WF4B (2 
WF4(2 , : , 1 , : ) =WF4B (2 , : , 

WF4(3, 1 , : , : )=WF4B(3, : , : ) ; 

WF4 (4,1,:,:) =WF4B (4 , : , : ) ; 

WF4(5, : , 1 , : ) =WF4B (3 , : , 

WF4(6 , : , 1 , : ) =WF4 (4,1, : , 

Ul  Initial  conditions 
for  L=l:Natoms 

for  m=2:Nmeas 
if  L>1 

WF4(3,L,l,m)=WF4(3,L-l,l,m-l)  . . . 

+delt . *HamEf fx3(3 , : ) *squeeze(WF4( : ,L-1 , 1 ,m-l) ) ; 
WF4(4,L,l,m)=WF4(4,L-l,l,m) . . . 

+delt . *HamEf fx3(4, : ) *squeeze(WF4( : ,L-1 , 1 ,m) ) ; 


end 
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if  RydInd(L, 1)==0 

WF4(7,L,l,m)=WF4(7,L,l,m-l)  +  delt . *  ... 

HamEf fx3(7, : ) ^squeeze (WF4( : ,L, 1 ,m-l) ) ; 

WF4(8,L, 1 ,m)=eitR/3 . *sum(squeeze (WF4 ( [4  6  7] ,L, 1 ,m) ) ) ; 

end 

end 

end 

m  Initial  conditions 
WF4(7, 1 , : , :)=WF4(7, : ,1,  :); 
m  Initial  conditions 
for  k=2:Natoms 

for  m=2:Nmeas 

WF4(6,l,k,m)=WF4(6,l,k-l,m) . . . 

+delt . *HamEffx3(6 , : ) ^squeeze (WF4( : , 1 ,k-l ,m) ) ; 
WF4(5,l,k,m)=WF4(5,l,k-l,m-l) . . . 

+delt . *HamEffx3(5, : ) ^squeeze (WF4( : , 1 ,k-l ,m-l) ) ; 


if  Rydlnd(l ,k)==0 

WF4(7, 1 ,k,m)=WF4(7, 1 ,k,m-l) . . . 

+delt . *HamEf fx3(7 , : ) *squeeze(WF4( : , 1 ,k,m-l) ) ; 

end 

if  Rydlnd(l,k)  ==  0 

WF4(8, 1 ,k,m)=eitR/3 . *sum(squeeze(WF4( [4  6  7] , 1 ,k,m) ) ) ; 

end 

end 

end 

m  Full  WF 
for  E=2:Natoms 

for  k  =  2:Natoms 
for  m=2:Nmeas 

WF4(2,L,k,m)=WF4(2,L-l,k-l,m)  . . . 

+delt . *HamEf fx3(2 , : ) *squeeze(WF4( : ,L_1 ,k-l ,m) ) ; 

WF4(6,L,k,m)=WF4(6,L,k-l,m) . . . 

+delt . *HamEf fx3(6 , : ) *squeeze(WF4( : , L,k-1 ,m) ) ; 


WF4(l,L,k,m)=WF4(l,L-l,k-l,m-l) . . . 
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+delt . *HamEf fx3(l , : ) *squeeze(WF4( : ,L-1 ,k-l ,m-l) ) ; 
WF4(5,L,k,m)=WF4(5,L,k-l,m-l) . . . 

+delt . *HamEf fx3(5 , : ) *squeeze(WF4( : , L,k-1 ,m-l) ) ; 


WF4(4,L,k,m)=WF4(4,L-l,k,m) . . . 

+delt . *HamEf fx3(4, : ) *squeeze(WF4( : ,  L-l ,k,m) ) ; 

WF4(3,L,k,m)=WF4(3,L-l,k,m-l)  . . . 

+delt . *HamEf fx3(3, : ) *squeeze(WF4( : , L-l ,k,m-l) ) ; 
if  RydInd(L,k)==0 

WF4 (7 , L , k , m) =WF4 (7 , L , k , ra-1 ) . . . 

+delt . *HamEf fx3(7, : ) ^squeeze (WF4( : ,L,k,m-l) ) ; 


end 

if  RydInd(L,k)  ==  0 

WF4(8,L,k,m)=eitR/3. *sum(squeeze(WF4( [4  6  7] ,L,k,m) ) ) ; 

end 


end 

end 


end 


1111111  Region  5  1111111111 


HI  Initial  conditions 

WF5 ( 1 , : , 1 , : ) =WF4 ( 1 , : , end , : ) ; 

WF5 (2 , : ,1, : ) =WF4 (5 , : .end, :) ; 

WF5(1, : , : , 1) =WF5 (1 , : , 1 , : ) ; 

WF5 (2 , : , : , 1) =WF5 (2 , : ,1,  :); 

WF5 ( 1 , 1 , : , : ) =WF5B ( 1 , : , : ) ; 
for  k=2:Nraeas 

for  m=2:Nmeas 

WF5(2,l,k,m)=WF5(2,l,k-l,m-l)+delt.*. . . 
HamEf f (2 , : ) ^squeeze (WF5 ( : , 1 ,k-l ,m-l) ) ; 

end 

end 

trap=squeeze (WF5 ( : , : , 1 ,  : ) )  ; 
for  k=0:Nmeas-2 

WF5diagk=zeros (2 , Natoms , (Nmeas-k) ) ; 
tmp=squeeze (WF5( : , : , 1 ,k+l) ) ; 
for  m=2 : (Nmeas-k) 

for  mm=l : round (c/vg) 
tmptmp=tmp ; 
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tmp(l , 2 : end)=tmptmp(l , 1 : end-1) +delt . * . . . 

HamEf f (1 , : ) *tmptmp( : , 1 : end-1) ; 
tmpC 2 , 2 : end)=tmptmp(2 , 2 : end)+delt . * . . . 

HamEf f (2 , : ) *tmptmp ( : , 2 : end) ; 

end 

WF5diagk(: , : ,m)=tmp; 

end 

for  L=2:Natoms 

WF5(1 ,L> : , : )=squeeze(WF5 (1 , L, : , : ) )+diag( squeeze (WF5diagk(l ,L, : ) ) ,k) ; 

WF5(2 ,L> : , : )=squeeze(WF5 (2,L, : , : ) )+diag( squeeze (WF5diagk( 2, L, : ) ) ,k) ; 

end 

end 

for  L=2:Natoms 

WF5(1  ,L,  :  ,  :  )=squeeze  (WF5(1 ,L,  :  ,  : )  )+squeeze  (WF5(1 ,L,  -diag(diag (squeeze (W 

end 

for  k=2:Nmeas 
tmptmp=tmp ; 
for  L=2:Natoms 

tmpCl ,  L,k: end)=squeeze(tmptmp(l ,L-1 , (k-1) : (end-1) ))+ . . . 

(delt . *HamEff (1 , : ) ^squeeze (tmptmp( : ,  L-l , (k-1) : end-1) ) ) . 5 ; 
tmp(2,L,k:end)=squeeze(tmptmp(2,L, (k-1) :end-l))+  . . . 

(delt*HamEf f (2, : ) ^squeeze (tmptmp( : ,L> (k-1) : end-1) ) ) . ’ ; 


end 

WF5 ( 1 , 2 : end , k , k : end) =tmp (1,2: end , k : end) ; 
WF5 ( 2 , 2 : end , k , k : end) =tmp (2 , 2 : end , k : end) ; 
tmp=WF5( : , : ,k, : ) ; 

end 

r///0  3  Photon  WF  or  equivalently  g3 
EEEmeas=squeeze (WF5 ( 1 , end ,:,:)); 
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